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NOTATION 


a  Duct  chord 

b(z)  Annular-airfoil  inner-surface  ordinate 

c(z)  Mean  line  ordinate  of  the  duct  section  measured  from  the  nose- 

tail  line 

1  2 

0^  [pfx^,  z)  “  P0]/yP  V  ,  Pressure  coefficient 

E(k)  Complete  elliptic  integral  of  the  second  kind 

h  (a/2  R^)  chord-diameter  ratio  of  the  duct 

K(k)  Complete  elliptic  integral  of  the  first  kind 

k  Modulus  of  the  elliptic  integrals 

z)  Local  pressure  on  the  annular  airfoil 
p^  Ambient  pressure  at  infinity 

q  Ring-source  strength 

Rd  Duct  radius 

(r,  <f),  z)  Cylindrical  coordinates 

s(z)  Half-thickness  ordinate  of  the  duct  section 

u(z)  Annular- airfoil  outer-surface  ordinate 

V  Free-stream  velocity 

wa  Axial  component  of  induced  velocity 

wq  Component  of  free-stream  velocity  in  direction  of  duct  axis 

Radial  component  of  induced  velocity 

w  Taylor  wake  fraction  at  the  duct 

xd 

x  Radial  coordinate  nondimensionalized  by  the  propeller  radius 

a  Angle  of  attack  of  a  duct  section 

y  R.ing-vortex  strength 

p  Mass  density 


IV 


Subscripts 

d  Duct  q  Ring  source 

p  Propeller  y  Ring  vortex 

Note:  Many  functions  are  defined  in  the  text. 
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ABSTRACT 

yk  computer  program  is  presented  that  calculates  the  annular 
airfoil  shape  from  a  given  pressure  distribution,  A  brief 
review  of  the  theory  of  the  invei'se  problem  of  the  annular  air¬ 
foil  is  also  presented.  The  distortion  of  the  duct  shape  by 
the  presence  of  an  axisymmetric  body  or  a  propeller  may  be 
taken  into  consideration.  Calculations  show  that  for  a  given 
pressure  distribution,  the  propeller  loading  and  location 
affect  the  duct  shape. 

ADMINISTRATIVE  INFORMATION 


This  work  was  covered  by  Subproject  S-R011  01  01,  Task  0401,  under 
the  Bureau  of  Ships  In-House  Independent  Research  Program. 

1.  INTRODUCTION 

Annular  airfoils  have  been  used  for  many  years1  as  shroud  rings 
around  marine  propellers.  These  ducts  have  been  used  for  the  acceleration 
of  the  velocity  at  the  propeller  (Kort  nozzles);  this  has  the  effect  of 
increasing  the  efficiency  when  the  unit  is  heavily  loaded.  Impetus  in  the 
last  few  years  has  been  given  for  their  use  to  increase  the  thrust  of  pro¬ 
pellers  during  takeoff  for  short  and  vertical  takeoff  aircraft.  Also,  in 
naval  architecture,  ducts  that  decelerate  the  flow  at  the  propeller  (pump- 
jets)  have  found  application  for  delaying  cavitation  on  the  propeller. 

A  number  of  experimental  and  theoretical  investigations  have  been 

conducted  on  the  annular  airfoil  and  the  ducted  propeller.  A  general 

2 

review  of  these  studies  has  been  made  by  Sacks  and  Burnell  so  only  the 

pertinent  and  more  recent  results  will  be  reviewed  here.  In  most  of  the 

theoretical  investigations,  a  distribution  of  ring  vortices  and  ring 

sources  is  used  which  lie  on  a  cylinder  of  diameter  representative  of  the 

3-7 

duct  and  of  length  equal  to  the  duct  length.  The  use  of  this  mathe¬ 
matical  model  implies  that  the  boundary  conditions  are  linearized  and  are 
satisfied  on  the  representative  cyliTider  and  not  on  the  duct  surface.  A 
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References  are  listed  on  page  46. 
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nonlinear  theory  for  the  annular  airfoil  has  been  presented  by  Chaplin 
where  the  annular  airfoil  is  represented  by  a  system  of  ring  vortices 
lying  on  the  surface  of  the  airfoil.  He  includes  both  the  static  and  the 
free- flight  cases. 

In  the  references  just  given,  the  direct  problem  of  the  annular 
airfoil  is  considered,  i.e.,  given  an  annular  airfoil  shape,  determine  the 
pressure  distribution  and  forces  on  the  duct.  Reference  9,  however,  pre¬ 
sents  the  theory  for  the  inverse  problem,  i.e.,  given  a  pressure  dist.ri- 
bution,  determine  the  annular  airfoil  shape. 

Both  the  direct  and  the  inverse  problems  require  the  solution  of  a 
singular  integral  equation,  the  first  for  the  ring  vortex  distribution, 
and  the  second  for  the  ring  source  distribution.  Another  approach  given 
in  References  10,  11,  12,  and  13  is  to  assume  the  ring  vortex  strength 

and,  if  the  effect  of  thickness  is  considered,  to  assume  the  thickness 
distribution.  The  disadvantages  of  this  procedure  are  that  the  circulation 
must  be  specified,  which  is  not  a  physical  property  of  the  annular  airfoil, 
and  it  is  not  possible  to  tell  a  priori  whether  the  pressure  distribution 
or  the  shape  will  be  satisfactory. 

The  usefulness  of  the  inverse  problem  is  to  delineate  annular  air¬ 
foil  shapes  which  will  operate  satisfactorily  for  a  given  flow  condition. 
This  is  to  say  that  in  the  presence  of  a  propeller  producing  a  given 
thrust,  a  duct  shape  can  be  determined  which  will  not  separate  in  air  or 
water,  nor  cavitate  in  water.  Both  separation  and  cavitation  are,  of 
course,  real  fluid  effects,  so  some  criteria  for  the  pressure  distribution 
must  be  specified  which  will  indicate  satisfactory  operation.  For  in¬ 
stance,  for  cavitation  a  minimum  pressure  coefficient  would  be  assumed, 
and  for  separation  a  maximum  rate  of  change  of  the  pressure  coefficient 
would  be  assumed. 

This  report  presents  a  computer  program  and  representative  type 
calculations  based  on  the  inverse  problem  presented  in  Reference  9.  The 
program  allows  the  inclusion  of  an  axisymmetric  perturbation  velocity  so 


There  is  a  difference  in  definition  of  the  direct  and  inverse  problem 
between  References  9  and  11. 
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that  the  duct  shape  can  be  determined  in  the  presence  of  a  propeller  or 
any  body  which  is  axisymmetric  to  the  duct.  Only  the  average  effect  of 
the  propeller  can  be  considered,  however,  as  including  the  effect  of  a 
finite  number  of  blades  would  imply  that  the  annular  airfoil  can  change 
shape  as  the  propeller  rotates. 

The  following  discussion  is  divided  into  three  main  sections.  The 
development  of  the  theory  is  reviewed  briefly,  then  the  computer  program 
is  described  and  presented,  and  finally,  computed  duct  shapes  are  pre¬ 
sented  for  comparison.  Calculations  are  also  included  which  show  how  the 
shape  of  the  duct  changes  with  variation  in  propeller  loading. 

2.  LINEARIZED  THEORY  OF  THE  ANNULAR  AIRFOIL 

Although  the  linearized  theory  of  the  annular  airfoil  has  been 
adequately  developed  in  the  reference  cited,  a  brief  sketch  of  the  develop¬ 
ment  will  be  repeated  here  for  the  sake  of  completeness. 

As  was  stated  in  the  Introduction,  the  method  of  singularities  is 
used  for  the  representation  of  the  flow  field  about  the  annular  airfoil. 

The  mathematical  model  used  is  a  distribution  of  ring. vortices  and  ring 
sources  lying  on  a  cylinder  having  a  length  equal  to  the  length  of  the 
annular  airfoil  and  a  diameter  representative  of  the  diameter  of  the 
annular  airfoil.  The  use  of  this  model  for  the  inverse  problem  necessitates 
a  number  of  assumptions: 

1.  The  real  fluid  is  inviscid  and  incompressible  and  no  separation 
occurs  on  the  annular  airfoil. 

2.  Body  forces  such  as  gravity  are  negligible. 

3.  The  freestream  flow  is  axisymmetric  and  steady,  but  an 
axis?mmetric  disturbing  velocity  may  exist.  An  implication  here  is  that 
the  static  condition  is  not  considered. 

4.  The  annular  airfoil  is  axisymmetric  and  has  finite  length. 

5.  The  distribution  of  ring  vortices  and  ring  sources  along  a 
cylinder  does,  indeed,  represent  the  annular  airfoil.  This  implies  that 
the  boundary  conditions  are  linearized. 
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The  coordinate  system  used  is  a  cylindrical  system  (r,  <}>,  z)  whose 
polar-axes  are  located  at  the  trailing  edge  and  z-axis  coincides  with  the 
centerline  of  the  annular  airfoil.  For  convenience,  the  axial  coordinate 
is  nondimensionalized  by  the  chord  a,  and  the  radial  coordinate  by  the 
representative  radius  R^,  of  the  annular  airfoil.  Figure  1  gives  the 
coordinate  system  and  Figure  2  is  a  delineation  of  the  system. 

The  pressure  distribution  on  the  annular  airfoil  is  assumed  known 
and  the  cross  section  of  the  airfoil  is  calculated  in  terms  of  a  camber 
distribution  c(z),  a  thickness  distribution  s(z),  and  an  angle  of  attack 
a,  as  follows: 


u^(z)  =  c'(z)  +  s'(z)  +  tan  a 
b^(z)  =  c'( z)  -  s^(z)  +  tan  a 


[2.1.1] 


where  u'(z)  is  the  slope  of  the  outer  surface  and  b"(z)  is  the  slope  of  the 
inner  surface  of  the  annular  airfoil. 

The  boundary  conditions  to  be  satisfied  are  that  the  normal  veloc¬ 
ity  must  be  zero  on  the  surface  of  the  annular  airfoil  and  the  Kutra  con- 

14 

dition  must  be  satisfied  at  the  trailing  edge.  In  linearize!  theory 
this  means  that  the  radial  velocity  at  the  representative  cylinder  must  be 
equal  to  the  slope  of  the  section,  or 


w 

— ^-(xd  +  0,  z)  =  -(1  -  wx  )  [c'(z)  +  s'(z).  +  tan  a]  [2.1.2] 

d 

w 

where  for  convenience  for  the  naval  architect  the  wake  (1  -•  w  )  =  ~  is 

Xd  V 

introduced.  At  the  trailing  edge  the  radia]  velocity  is  zero,  or 


w 

-y?-  Od  ±  0,  0)  *  0 

In  the  symbol  +,  the  +  sign  denotes  the  outer  surface  of  the  annular  air¬ 
foil  and  the  -  sign  denotes  the  inner  surface. 
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2.2  DERIVATION  OF  THE  EQUATION  FOR  THE 
THICKNESS  DISTRIBUTION 

To  calculate  the  thickness  distribution,  camber  distribution,  and 
the  angle  of  attack  of  the  annular  airfoil  from  the  pressure  distribution, 
it  is  necessary  to  obtain  the  strength  of  the  ring  vortices  and  ring 
sources.  The  ring  vortex  and  ring  source  strengths  must  be  of  sufficient 
magnitude  that  they  induce  radial  velocities  that  satisfy  the  given 
boundary  conditions.  By  substituting  the  equations  for  the  radial  veloc¬ 
ity  induced  by  the  ring  vortices  and  sources  into  Equation  [-2.1.1],  an 
expression  for  the  slope  of  the  thickness  distribution  is  obtained  in 
terms  of  known  quantities. 

The  nondimensional  elementary  circulation  of  the  ring  vortices  is 
represented  by  yCz)  and  the  strength  of  the  ring  sources  by  q(z).  Since 
the  velocities  are  linear,  they  are  additive;  thus  the  radial  velocities 
along  the  annular  airfoil  are  found  to  be 


w 

~~  Ud  ±  0,  z)  = 


V  Cxd  1  °>  Z) 


V  Z 


(0  <  z  <  1) 


IP 

[2.2.1] 


[wr  (xd,  z)]^  =  the  radial  velocity  induced  on  the  annular 

airfoil  by  the  ring- vortex  system, 

[wr  (xd  +  0,  z)]  =  the  radial  velocity  induced  on  the  annular 

airfoil  by  the  ring-source  system,  and 


[wr  (xd,  z)]^  =  the  radial  velocity  induced  on  the  annular 

airfoil  by  a  propeller  or  any  other  singu¬ 
larities  in  the  flow. 

The  expressions  for  these  velocity  components  have  been  derived  in 
Reference  14  and  are  as  follows: 


fK(k)  -  E(k)]  -  2E(k)}dz' 

[2.2.2] 
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'  W  "I  ,  . 

~  (xd  t  0,  z)Jq.=  r  q(z-)k  [K(k]  -  B(k)]  dz'  +  -±-  q(z) 


12.2.3] 


where  h  =  k2  =  — 


—x - « - ,  and  K(k)  and  E(k)  are  complete 

li  (z  -  Z'y  +  1 


elliptic  integrals  of  the  first  and  second  kind,  respectively.  The  symbol 
J  denotes  a  Cauchy  principal  value  integral. 

By  substituting  these  expressions  in  Equations  [2.1.2]  and  [2.2.1], 
the  following  equation  is  obtained. 


x 

-h  f  k  {4h20-z")2  tK(k)  -  E (k) ]  -  2E(k)}  dz- 

0 

+  4 :  f1  q(z')k  [K(k)  -  E(k)]  dz"  +  \  q(z)  = 

17  0 

"  (1_Wxdj  -  S'Q  +  tan  -  y  ^xd>  z)  p 


[2.2.4] 


Since  tnc  integrals  occurring  in  this  equation  have  only  one  sign  and  since 
the  radial  velocity  induced  by  the  propeller  on  the  annular  airfoil  does 
not  change  sign  from  one  side  of  the  foil  to  the  other,  it  may  be  concluded 
that  the  +  signs  go  with  the  +  signs  and  hence 


q(z)  =  -  2(1  - 


(l  ■  \) 


s'(z) 


[2.2.5] 


This  equation  gives  the  relationship  between  the  ring  source 
strength  and  the  slope  of  the  thickness  distribution.  The  next  step  is  to 
find  the  total  velocity  tangent  to  the  airfoil  in  terms  of  the  ring  vortex 
strength,  ring  source  strength,  and  the  strength  of  any  other  singularity 
present  in  the  flow  such  as  a  propeller.  Since  the  annular  airfoil  is 
theoretically  replaced  by  a  cylinder,  the  velocity  in  question  is  just  the 
axial  velocity.  The  pressure  distribution  is  related  to  the  induced  ve¬ 
locities  by  means  of  the  linearized  Bernoulli  equation,  i.e.. 


„  p(xd>  Z)  -  ,  wa 

P  ‘  1  ,,2  “V 

2  P  V 


y  (Xd>  Z) 


[2.2.6] 
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where  pCx^,  z)  =  the  local  pressure  on  the  annular  airfoil, 
pm  =  ambient  pressure  at  infinity, 

Wa(Xd'  =  axaa*  component  of  induced  velocity  on  the  annular 

airfoil,  and 

p  =  the  mass  density. 


Within  linearized  theory  all  velocities  are  additive;  thus  the 


total  axial  velocity  is  given  by 

’"a1  rvz)l  .  K(z)]  [wa(z)' 

V~Jtotal  ”[■  V  JT\  V  jq  +  [  V  jp 


w 

o 

V 


[2.2.7] 


=  the  axial  velocity  induced  on  the  annular  airfoil 
by  the  ring-vortex  system, 

-  the  axial  velocity  induced  on  the  annular  airfoil 
by  the  ring-source  system, 

=  the  axial  velocity  induced  on  the  annular  airfoil 
by  a  propeller  or  any  other  singularity  in  the  flow, 
and 

=  speed  of  advance  of  the  duct. 


From  Equation  [2.2.2]  the  perturbation  velocity  distribution  w/  is  obtained 

d 

as 


It  is,  of  course,  not  possible  to  determine  the  section  shape  in 
the  presence  of  a  propeller  with  a  finite  number  of  blades  since  the 
pressure  distribution  is  essentially  time-dependent;  however,  it  is  pos¬ 
sible  to  consider  the  average  effect  of  the  propeller  which  corresponds  to 
an  infinitely  bladed  propeller.  Under  this  restriction  the  duct  circu¬ 
lation  distribution,  defined  as  y(z),  is  a  function  of  z  only  and  justifies 
the  use  of  the  following  relationship  from  Reference  14  for  the  axial 
velocity  induced  on  the  duct  by  the  ring  vortex  distribution  and  ring 


8 


■StSisvc— 


source  distribution,  respectively: 
w  (x,  «•  0,  z)1  ,  1 

—  v' - L  =  §r/  y(z')  k  CK(k)  -  E (k) ]  dz'  +  f  y(z)  [2.2.9] 

L  v  JY  0  z 


.i(WkEmh. 

L  V  Jq  it  J  (z-z')  ^ 


[2.2.10] 


or  in  terms  of  Equation  [2.2.5] 


\(v  z) 

M 


/1-w  \  1 

® k  E(k)d2' 


If  Equations  [2.2.9]  and  [2.2.10]  are  substituted  into  Equation 


[2.2.8], 


wa(z)  1  h  1  1 

-y-J.  =icP  -it'  k  tK(k)  -  E<k)i dz' + 


t1  Wxd)  A  s'(z')  k  E(1Q  .  ,  wa(z) 

»  J  (*-*')  "  v 


V  Ip 


[2.2.11] 


A  mean  velocity  is  defined  as 


wa00 


i  Kc*)i  rw.wi  i 

=  T \  v  J  +  +  L^T-J-J  =  T [V  +  cp-]  [2'2a2] 


V  mean  2  I  V 


and  substitution  of  Equation  [2.2.11]  into  this  equation  gives 


w  (z)l  ,  } 

~V  jmean  S  /  ^  k 


[K(k)  -  E (k) ]  dz ' 


[2.2.13] 


/1-w  \  1 

l  Xd /  £  s'(z')  k  E 
■  it  J  (z  -  z' 


lULgw  dz'  ♦ 

Z  -  z')  v 


V  IP 


A  velocity  difference  is  also  defined  as 


v  (z) 
ak  _ 

V  diff 


wa(z)|  |waCz) 


*  v  j-r  2  [ 


C  A  -  C 

p+  p- 


[2.2.14] 


and  substitution  into  Equation  [2.2.11]  gives 
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~1 


4 

'b 


waC^l 


V 


diff 


Y(z) 


[2.2.15] 


Since 

w  (z) 

and 

’V2) 

j 

V 

mean 

L  V 

diff  are  known  from  t^ie  pressure  distribution, 
substitution  of  Equation  [2.2.15]  into  [2.2.13]  gives  a  singular  integral 
equation  for  the  slope  of  the  thickness  distribution  in  terms  of  known 
quantities : 


A 

f 


S'(z)  k  C(k)  _  -  7T  l 


w  (z) 


mean 


Vz) 


pj 


V2) 


[2.2.16] 


di££  k  [K(k)  -  E(k) ]  dz' 


2.3  REDUCTION  OF  THE  THICKNESS  DISTRIBUTION 
EQUATION  FOR  SOLUTION  ON  THE  IBM  7090 

To  facilitate  solving  the  expression  for  the  slope  of  the  thickness 
distribution.  Equation  [2.2.16]  is  reduced  to  a  Fredholm  equation  of  the 
second  kind  by  using  a  method  given  by  Muskhelishvili. 15  The  resulting 
equation  contains  singularities  and  is  handled  by  a  method  discussed  in 
Reference  16. 

For  convenience,  Equation  [2.2.16]  is  rewritten  as 


£  z(z ' z' rt'-Kh  dz' =  HCz) 


[2.3.1] 


where  g(z  -  z')  =  k  E(k) 
and  H(z)  =  - 


h 

i 

r 

'V2')" 

T\)  i 

i  - 

i 

V 

dif£  k  [*00  -  E(k) ]  dz'  - 


\(Z)' 

w  (z) 
av 

V 

mean 

V 

P> 

r-\) 

s'fz'] 

Now  the  term  g(0)  -  -~j  is  added  to  and  subtracted  from  the  integrand  of 

Equation  [2.3.1]  giving 
1 

SCO)  j  dz'  +  C1  [g(z  -  Z-)  -  8(0))  dz '  =  H(z) 
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^22 


or  since  g(0)  =  1 


a 

f  dz'  =  H(z)  -  Z1  [g(z  -  z')  -  1]  dz'  =  H(z)  [2.3.2] 

o_  0 

This  equation  is  now  in  the  form  of  a  Cauchy  integral  equation  because  the 
right  side  is  free  from  singularities.  As  found  in  Reference  9,  this 
equation  has  a  unique  inverse  given  by 


s'(z)  = 


v.  vz  (1 


i — fi  f  fl(z-) 

5~ T)  I”  J  -  z) 


dz'  +  2  Z1  s'(z')  dz' 


Since  the  airfoil  section  must  be  closed,  /  s'(z')  dz'  =  0  and  by  using 

0 

equation- [2.3.2]  a  Fredholm  equation  of  the  second  kind  is  obtained 


s'(z)  = 


where 


z  (1  -  z) 


f(z)  +— -Z1  G(z,-  z')  s'(z')  dz' 
ir  0 


[2.3.3] 


f(2)  -  f  d2' 


G(z,  z') 


i  /  14"  (l  -  z" 
~  2  J  (.*"  ~  z) 


)  [1  -  g(z'  -  z")] 


Lz  -  z 


Equation  [2.3.3]  cannot  be  handled  in  its  present  form  since  there  exists 
a  square  root  singularity  at  z  =  0  and  z  =  1.  To  remove  the  singularities 
and,  as  will  be  seen  later,  to  facilitate  evaluation  of  f(z)  and  G(z,  z') 
the  following  change  of  variable  is  used: 


Since 


z  =  Y  (1  +  cos  0) 


,  2  s ' [(1/2)  (1  +  cos  6)] 

s  W  sliTe 


the  following  equation  is  obtained  from  Equation  [2.3.3]: 


s*'(e)  =  f*  (0)  +  f1  G* (9,9')  s* ^ (0 ')  d0 ' 

0 


[2.3.4] 
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-Vw 


where 


f\e)  =  .  I _  /  H  (6'*)  sin2  8 "  dQ'* 

2tt2  J  (cos  8'  ~  cos  0) 


G*(G,  e-)  =JL  r  .II.-  fc  E(k)l  Sin2  6~  d0~ 
2  /  fcos  ( nr*. c  a -* ^ 


tt2  j  (cos  Q"  -  cos  6)  (cos  0""  -  cos  0") 


k2  = 


1  ,2  ,  > 

X"  (cos  6~  -  cos  Q')  +  1 

The  symbol  *  is  used  to  simplify  notation;  for  example, 

f(z)  =  f[(l/2)  (1  +  cos  0)]  =  f*(e) 

The  integrals  f  (0)  and  G  (0,  0-)  are  Cauchy  principal  value  inte- 

grals  and  to  evaluate,  part  of  the  integrand  of  each  is  expanded  in  a 
Fourier  cosine  series,  i.e., 


where 


H  (0')  sin2  d'  =  a  +  I 


0  n  —  i  «v 


a0  -7/*  H*ceo  sin2  e-  de' 


a  cos  0' 


then 


an  h  (0')  sin2  0'  cos  n0'  d0' 


f*(8)  =  4-  £  a 

2^  n_2  n  sm  0 


[2.3.5] 


where 


[1-k  E(k)]  sin2  e"  u  ,  „ 
con  e"  "cos  6'  =  b0(9  J  +  l  b  (0')  cos  n0- 

n=i 


Ve'>  -  7  / 


[1-k  E(k) 1  sin2  e 


cos  d"  ~  cos  0 


7—  d0- 


then 


,  r«^  -  2_  f  [1-k  E(k)3  sin2  Q" 

;  7T  J  Tos  Q"  -  cos  0 " . .  cos  n6~  d0"' 
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G*(0,  0")  =  -  Z  b  (60  S1? 
v  ’  7T  .  nv  '  sm  0 

n=l 


[2.3.6] 


The  kernel  G  (0,  0")  is  now  of  the  degenerate  type  and  Equation  [2.3.4] 
can  be  solved  by  the  method  applicable  to  this  type  of  kernel.  Details  of 
the  method  are  given  in  Reference  17.  If  this  procedure  is  followed,  the 
following  equation  for  the  slope  of  the  thickness  distribution  is  obtained. 


1  .  sin  26 


"  fn-,  J?"  rn\  l  a  .  i  .  Sin  .  1  A  Sin  n0 

S  (0)  «  f  (0)  +  —  A.  +  A0 - 7 T-  +  ...  +  A  - 7 — — 

K  v  J  ir  1  ii  2  sm  0  n  n  sm  0 


where 


•k 

f  (0)  is  given  by  Equation  [2.3.5]. 

And  Afi  is  given  by  the  following  set  of  simultaneous  equations: 

A1  C1  ”  Cll)  '  A2  C12  ~  ~  An  Cln  =  dl 

-  A1  c21  +  A2(l  -  c22)  -  . “  An  c2n  =  d2 


[2.3.7] 


[2.3.8] 


where 


-  K1  Cnl  -  A2  cn2  '  .  +  An  f1 


)  =  “n 


c. .  =—f'  b.(6')  51”  do-  (i,1  =  1,2.3  ...  n) 

v  0  iK  J  sm  0  v  J 

d.  =  /*  b.  (0")  f*(6')  d0'  (i  =  1,2,3  ...  n) 

1  0  1 

Expression  [2.3,7]  is  the  expression  solved  on  the  computer  with  a 

finite  sum  replacing  the  infinite  sum  in  Equation  [2.3.5].  The  c^’s  and 

d^'s  may  be  evaluated  as  given,  so  the  AJs  follow  immediately  from 

[2,3.8].  However,  tc.1  evaluate  the  a  's  in  Equation  [2.3.5],  the  singularity 

*  n 

must  be  removed  in  H  (0').  This  is  accomplished  by  using  the  change  of 
3 

variables  t  =  Q"  -  Q'  where  Q"  is  the  variable  of  integration  for 
* 

H  (0').  The  thickness  distribution  now  follows  immediately  by  integrating 
Equation  [2.3.7]  from  0  to  0. 
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2.4  DERIVATION  OF  THE  EQUATION  FOR  THE  CAMBER 
DISTRIBUTION  AND  ANGLE  OF  ATTACK 

To  calculate  the  camber  distribution  and  the  angle  of  attack,  it  is 
only  necessary  to  integrate  Equation  [2.2.4]  after  substituting  for  the 
known  quantities.  Equation  [2,2.4]  is  rewritten  and  substitution  for  the 
source  strength  q  and  vortex  strength  y  is  made  from  Equations  [2.2.5]  and 
[2.2.15],  respectively. 


c"(z)  +  tan  <x  = 


.2  l  fw  ( z ') 

Tjihr  J  [nr-jdiff (z  - z^  k  [K(k)  - E(k)]dz' 

V  Xdi  0 


+ _ 1  f  [”a(l°]  k  E(k) 

2ir|l  -  wx  j  J  |_  v  Jdiff  (z  -  z ') 

-  swi"*'-  Trhn  L^(v2)J, 


[2.4.1] 


If  Equation  [2.4.1]  is  integrated  with  respect  to  z  from  0  to  1, 

the  tangent  of  the  angle  of  attack  is  obtained  since  f1  c'(z)  dz  =  0. 

0 

Equation  [2.4.1]  is  the  desired  expression  for  the  slope  of.  the  camber 
distribution. 


2.5  SIMPLIFICATION  OF  THE  CAMBER  DISTRIBUTION 
EQUATION  FOR  SOLUTION  ON  THE  IBM  7090 

To  facilitate  solving  the  expression  for  the  slope  of  the  camber 
distribution,  Equation  [2.4.1],  a  change  of  variable  is  made.  The  resulting 
equation  contains  various  singularities  which  are  handled  by  a  method 
suggested  in  Reference  15.. 

With  z  =  y  (1  +  cos  0]  Equation  [2.4.1]  has  the  form 


, -2c  (8) 

sin  0 


+  tan  a  = 


Jl _ f 

f1  -  W*.U 


r 


d/  0 


diff 


(cos  0  -  cos  0") 


k[K(k)  -  E(k)]  sin  0'  d6' 


t  K*ce) 


2ir/l  -  w  \ 

l  xd) 


k  E(k)  sin  6"  de' 
V  diff  cos  0  -  cos  0" 


a*., 


_  A  fv  s*\q')  k[K(k)  -  E(k)]  d0' 
v  0 


Wr  (xd>0) 


[2.5.1] 


(1  -  wXd)  L  v  Jp 

where  the  notation  *  and  the  form  of  k  are  given  in  Equation  [2.3.4]. 

Because  of  the  presence  of  the  term  cos  0  -  cos  Q'  the  first  inte¬ 
gral  in  Equation  [2,5.1]  has  no  singularities  and  may,  therefore,  be 
evaluated  immediately  by  use  of  the  computer. 

The  second  integral  in  Equation  [2.5.1]  is  a  Cauchy  principal  value 
integral  and  is  first  simplified  by  a  technique  given  in  Reference  15;  the 


term 


w  *(9) 


_  sin  0 

V  diff  cos  0  -  cos  Q' 


is  added  to  and  subtracted  from  the  integrand: 


2ir/l  -  w 


-  k*(e) 

f  L  V  Jdiff 


kE(k)  sin  6 


w.  (9) 


-V— Jdiff  sin  9 


cos  0  -  cos  0' 


d  v  0 


r  *  IT 

wa  (0)  .  a  f  de' 

+  V  diff  sin  9  J  cos  0  -  co,'  0" 


[2.5.2] 


The  second  integral  in  this  expression  is  zero  so  only  the  first  integral 
needs  to  be  considered. 

The  behavior  of  the  integrand  is  investigated  as  0'  -*■  0  by 
L' Hospital’s  rule  and  the  following  limit  is  obtained: 


/  [Vw 

•’  de  IL  V  Jdiff 

sin  0 


sin  8 


[2.5.3] 


Now  with  this  knowledge  of  the  behavior  of  the  integrand  of  expression 
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[2.5,2],  the  second  integral  in  Equation  [2.5.1]  is  evaluated  by  inte¬ 
grating  [2.5.2]  on  the  computer  and  using  [2.5.3]  evaluated  numerically 
for  the  value  of  the  integrand  at  6"  =  0. 

The  third 'integral  in  Equation  [2.5.1]  is  simplified  as  follows. 
First  this  equation  is  integrated  by  parts  with 

u  =  k  [K(k)  -  E(k)] 


dv  =  s  '(0")  d9" 

* 

Then  v  =  s  (S')  and  after  differentiation  and  considerable  simplification 
(see  page  107  of  Reference  18) , 

du  =  -  k3  j-  (cos  0  -  cos  0 ")  sin  0"  •  K(k)  -  E(k)  +  E(k)  \  do' 

it  it  1  —  k  / 

Now  by  using  the  fact  that  s  (0")  L^_n  =  s  (6 ")  1  „ ^  =  0 

D  — U  0  =TT 

the  following  equation  is  obtained 

2 

s*'  (O')  k[K(k)  -  E(k)  ]  d@"  =x"/71  s*(0')  k'3  (cos  0  -  cos  0") 


sin  9'  [K(k)  -  E(k) ]  d0' 


J  co 


*  (8  ^)  k3  sin  9'  E(k) 
cos  0  -  cos  8' 


[2.5.4] 


The  first  integral  on  the  right-hand  side  of  Equation  [2.5.4]  can  be 
evaluated  on  the  computer  immediately;  the  second,  however,  must  be 

simplified  as  in  the  second  integral  of  Equation  [2.5.1].  Thus  the  term 

* 

s  (6)  sin  6 
cos  0  -  cos  0" 

is  added  to  and  subtracted  from  the  integrand  obtaining 


/ ^ 


(0*)  k3  sin  0^  E(k)  -  s  (8)  sin  6 
cos  0  -  cos  0" 


since,  as  before, 


cos  0  -  cos  0' 


d0'  =  0 


16 


0 


As  in  the  integral  [2.5.2]  the  behavior  of  the  integrand  as  0'  ->■  0 
is  investigated  and  s  (0)  +  s  (0)  cot  0  is  obtained  as  the  value  of  the 
limit.  Thus  the  complete  expression  [2.5.4]  may  now  be  evaluated  numer¬ 
ically. 

All  three  integrals  of  Equation  [2.5.1]  are  now  in  a  form  suitable 
for  numerical  evaluation,  so  tan  a  may  be  found  by  integrating  [2.5.1] 

IT  *' 

from  0  to  it  and  recalling  that  /  c  (0  )d0"  ■-  0.  The  slope  of  the  cam- 

0 

ber  distribution  then  follows  immediately,  and  integration  from  0  to  0 
gives  the  camber  distribution. 


3.  COMPUTER  PROGRAM 

The  calculations  based  on  the  theory  presented  in  Section  2  have 
been  programmed  for  the  IBM  7050  high-speed  computer.  Input  consists  of  a 
pressure  distribution  on  the  inner  and  outer  surfaces  of  the  annular  air¬ 
foil,  the  chord-diameter  ratio  of  the  annular  airfoil,  and  a  wake 
fraction.  It  is  also  possible  to  make  calculations  with  a  propeller  or 
any  axisymmetric  body  located  in  the  annular  airfoil;  however,  axial  and 
radial  components  of  the  induced  velocity  must  be  included. 

The  output  consists  of  the  thickness  distribution,  the  camber 
distribution,  and  the  angle  of  attack  of  the  annular  airfoil.  It  takes 
approximately  8  minutes  on  the  IBM  7090  high-speed  computer  to  obtain 
these  results. 

The  input-output  format  and  the  FORTRAN  listing  of  the  computer 
program  are  discussed  in  the  following  sections. 

3.1  INPUT  FORMAT 

The  first  input  card  contains  an  integer  which  represents  the  number 
of  cases  to  be  run.  This  number  is  entered  in  Columns  1-4  under  an  14 
format.  The  only  limit  to  the  number  of  cases  .to  be  run  is  a  consideration 
of  computer  time. 

The  second  card  is  a  problem  identification  card.  A  one  (1)  must 
appear  in  Column  1  and  any  alphanumeric  characters  may  appear  in  Columns  2 
through  72.  Even  if  no  identification  is  desired,  the  card  must  be  in¬ 
cluded  with  a  one  (1)  in  Column  1. 
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The  third  card  contains  10  parameters  in  the  format  2F10.6,8I4. 

They  are  as  follows: 

1.  The  chord- diameter  ratio. 

2.  The  wake  fraction. 

3.  The  maximum  number  of  Fourier  coefficients  to  be  allowed  in  the 

* 

evaluation  of  the  kernel  function  G  (0,8"),  Equation  [2.3.4].  The  user 
may  specify  any  number  not  exceeding  40.  The  program  calculates  Fourier 
coefficients  until  either  a  convergence  criteria  of  10”^  has  been  met  or 
the  number  of  coeff.‘  cients  equals  the  number  specified  in  this  parameter. 
The  program  normally  uses  between  20  and  25  terms  so  the  user  may  specify 
40  since  only  as  many  terms  are  calculated  as  are  needed. 

4.  The  number  of  ordinates  used  in  the  integral  evaluations.  The 
integrations  are  performed  using  Simpson's  rule  so  this  number  must  be 
odd.  The  maximum  allowed  by  the  dimension  of  the  program  is  101.  Nor¬ 
mally,  a  value  of  51  is  more  than  adequate  and  if  time  is  a  factor,  25 
would  probably  suffice. 

5.  The  number  of  input  pressure  data  points.  This  number  may  not 
exceed  37.  If  possible,  to  facilitate  interpolation,  the  user  should 
specify  closer  spaced  data  in  regions  where  the  pressure  curves  have  steep 
slopes. 

6.  The  number  of  points  to  be  used  for  the  harmonic  analysis  of 

★ 

the  Cauchy  principle  value  integral  f  (0),  Equation  [2.3.4].  The  maximum 
allowed  by  the  dimension  of  the  program  is  200  and  this  is  also  the 
suggested  value. 

7.  The  number  of  harmonics  used  in  the  harmonic  analysis  described 
in  parameter  6.  This  number  must  be  strictly  less  than  one-half  the  num¬ 
ber  specified  in  6.  Normally,  a  value  of  45  is  more  than  adequate. 

*  '** 

8.  The  maximum  number  of  ordinates  used  in  integrating  s  (0)  and 
c  (0)  to  obtain  the  thickness  and  camber  distributions.  The  number  must 
be  greater  than  or  equal  to  51  and  must  be  odd.  A  value  of  61  has  proven 
to  be  satisfactory.  There  is  no  upper  bound  on  this  parameter. 

9.  The  parameter  signifying  the  presence  of  a  propeller  or  any 
axisymmetric  body  in  the  annular  airfoil.  The  parameter  must  have  the 
values : 

1  if  a  propeller  or  body  is  present, 

2  if  a  propeller  or  body  is  not  present.  - 
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10.  The  parameter  controlling  the  punching  of  one-half  the  thick¬ 
ness  and  camber  distributions  on  IBM  cards  by  the  machine.  This  parameter 
must  have  the  values: 

1  if  no  data  is  to  be  punched, 

2  if  data  is  to  be  punched. 

The  format  of  the  punched  data  is  discussed  in  the  following  para¬ 
graphs.  A  discussion  of  the  effects  of  varying  parameters  3  and  7  may  be 
found  in  Reference  16. 

The  fourth  card  is  the  first  card  containing  pressure  data.  The 
value  of  parameter  5  on  the  third  card  determines  the  number  of  these 
cards.  Each  card  contains  the  abscissa  of  the  pressure  distribution  and 
four  pressure  data  terms  in  a  format  of  5F14.8.  The  fourth  card  contains 
the  abscissa  x  =  0.0  and  the  appropriate  pressure  data  at  the  leading  edge 
of  the  annular  airfoil.  The  fifth  card  contains  the  next  abscissa  and  its 
corresponding  pressure  data.  The  abscissa  value  may  be  arbitrarily  spaced 
in  the  interval  [0,1]  as  long  as  the  first  one  is  0.0  and  they  are  strictly 
increasing.  The  second  term  on  each  card  is  the  pressure  on  the  outer 
surface  of  the  annular  airfoil.  The  third  is  the  pressure  on  the  inner 
surface  of  the  annular  airfoil. 

If  no  propeller  is  present,  i.e.,  a  two  (2)  has  been  given  as 
parameter  9  on  the  third  card,  the  fourth  and  fifth  terms  on  each  card  are 
omitted;  however,  if  a  propeller  is  present,  i.e.,  a  one  (1)  has  been  given 
as  parameter  9  on  the  third  card,  then  the  fourth  and  fifth  terms  must  be 
specified.  The  fourth  term  is  the  axial  induced  velocity  and  the  fifth  is 
the  radial  induced  velocity.  An  .example  showing  the  input  data  for  a 
ducted  propeller  is  shown 'in  Appendix  A. 

All  cards  after  the  first  card  must  be  repeated  for  each  case  even 
though  some  of  the  data  may  be  the  same. 

3.2  OUTPUT  FORMAT 

The  f-'rst  page  of  output  contains  the  input  data.  The  second  page 
contains  the  angle  of  attack  of  the  annular  airfoil  in  degrees  and  a  table 
consisting  of  one-half  the  thickness  and  of  the  camber  distribution  from  the 
leading  edge  of  the  foil,  0.0,  to  the  trailing  edge,  1.0,  in  increments  of 
0.05.  The  output  obtained  from  the  input  data  of  the  ducted  propeller 
given  in  Section  3.1  is  also  shown  in  Appendix  A. 
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If  parameter  10  on  the  third  card  is  a  two  (2),  data  will  be 
punched  on  IBM  cards  by  the  program.  The  first  three  cards  contain  a 
total  of  21  numbers,  ranging  from  0.0  to  l.n  in  increments  of  0.05,  which 
are  the  abscissa  values  of  one-half  the  thickness  distribution.  The  next 
three  cards  contain  the  values  cf  one-half  the  thickness  distribution  at 
the  abscissa  values  given  the  first  three  cards.  The  next  three  cards, 
which  are  the  abscissa  values  of  the  camber  distribution,  are  the  same  as 
the  first  three  cards.  The  last  three  cards  contain  the  values  of  the 
camber  distribution  at  the  abscissa  values  given  on  the  preceding  three 
cards.  All  this  output  is  punched  in  a  format  of  9F8.5. 

3.3  FORTRAN  LISTING 

The  FORTRAN  listing  of  the  computer  program  is  given  in  Appendix  B. 
In  addition  to  the  subroutines  furnished  automatically  by  the  Bel-  Monitor 
System  on  the  IBM  7090,  the  binary  coding  of  the  following  subroutines 
available  from  SHARE  must  be  added  to  the  FORTRAN  listing:  BE-ELIP, 

AMGMHA,  AMMATI,  LACBRT,  AQALLI,  and  VG-AS  +  C.  The  program  takes  about 
eight  minutes  to  run  under  the  Bell  Monitor  System. 

4.  RESULTS  OF  CALCULATIONS 

4.1  ANNULAR  AIRFOIL 

A  number  of  calculations  were  made  to  compare  calculated  duct  shapes 
with  actual  shapes.  Two  annular  airfoils  from  Reference  7  were  con¬ 
sidered,  Duct  II  and  the  BTZ  duct,  and  the  pressure  distributions  for 
these  ducts  are  given  in  Tables  1  and  2,  respectively.  There  are  two 
pressure  distributions  given  for  each  duct;  one  is  the  experimental  dis¬ 
tribution  and  the  other,  the  linearized  theoretical  distribution  as  calcu¬ 
lated  by  the  method  of  Reference  7.  Tables  3  and  4  give  the  tabulated 
data  for  the  thickness,  camber,  and  angle  of  attack  as  calculated  and  the 
actual  shape  for  each  of  the  ducts.  Figures  3  and  4  show  a  comparison  of 
the  section  shapes  for  Duct  II  and  the  BTZ  duct,  respectively.  The 
ordinates  have  been  expanded  to  accentuate  the  differences. 

A  comparison  of  the  results  show  that  for  either  the  linearized  or 
experimental  pressure  distribution  the  calculated  thickness  is  a  few 
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TABLE  1 


Pressure  Distribution  for  Duct  II 


1  -  z 

Experimental 

Linearized  Theory 

°P  + 

C 

P- 

V 

C 

P- 

Leading  Edge 

1.000 

1.000 

1.000 

-1.000 

.  0.0019 

0.400 

o.4oo 

0.563 

-0.676 

0.0076 

0.151 

0.151 

0.117 

-0.247 

0.0170 

0.030 

0.030 

-0.077 

-0.062 

0.0302 

-0.08l 

0.058 

-0.195 

0.051 

0.0469 

-0.l60 

0.088 

-0.272 

0.121 

0.0670 

-0.199 

0.101 

-0.319 

0.162 

0.0904 

-0.217 

0.110 

-0.343 

0.179 

0.1170 

-0.234 

0.112 

-0.350 

0.178 

0.1464 

-0.250 

0.109 

-0.345 

0.166 

0.1786 

-0.275 

0.101 

-0.334 

0.148 

0.2132 

-0.305 

0.100 

-0.323 

0.129 

0.2500 

-0.321 

0:107 

-0.313 

0.114 

0.2887 

-0.331 

0.119 

-0.309 

0.107 

0.3290 

-0.338 

0.131 

-0.311 

0.108 

0.3706 

-0.339 

0.137 

-0.320 

0.118 

0.4132 

-0.345 

0.135 

-0.334 

0.134 

0.4564 

-0.351 

0.145 

-0.351 

0.155 

0.5000 

-0.355 

0.146 

-0.367 

0.177 

0.5436 

-0.350 

0.152 

-0.381 

0.198 

O.5868 

-0.360 

0.159 

-0.388 

0.215 

0.6294 

-0.369 

0.165 

-0.387 

0.227 

0.6710 

-0.359 

0.161 

-0.376 

0,233 

0.7113 

-0.345 

0.175 

-0.354 

0.233 

0.7500 

-0.330 

0.171 

-0.321 

0.230 

0.7868 

-0.308 

0.179 

-0.279 

0.224 

0.8214 

-0.270 

0.191 

-0.229 

0.218 

0.8536 

-0.235 

0.192 

-0.174 

0.213 

0.8830 

-0.169 

0.192 

-0.116 

0,209 

0.9096 

-0.100 

0.192 

-0.058 

0.209 

0.9330 

-0.032 

0.190 

-0.002 

0.211 

0.9532 

0.032 

0.190 

0.051 

0.215 

O.9698 

0.085 

0.192 

0.098 

0.220 

0.9830 

0.140 

0.195 

0.139 

0.224 

0.9924 

0.160 

0.197 

0.170 

0.224 

0.9981 

0.184 

0.199 

0.192 

0.218 

1.0000 

0.200 

0.200 

1.000 

1.000 
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TABLE  2 


Pressure  Distribution  for  the  BTZ  Duct 


1  -  z 


Experimental 


Linearized  Theory 


Leading  Edge 
0.0019 
0.0076 
0.0170 
0.0302 
0.0469 
0.0670 
0.0904 
0.1170 
0.1464 
0.1786 
0.2132 
0.2500 
0.2887 
0.3290 
0.3706 
0.4132 
0.4564 
0.5000 
0.5436 
0.5868 
0.6294 
0.6710 
0.7113 
0.7500 
0.7868 
0.8214 
0.8536 
0.8830 
0.9096 

0.9330 

0.9532 

0.9698 

0.9830 

0.9924 

0.9981 

1.0000 


1.000 

0.700 

0.120 


0.037 

0.085 

0.100 

0.110 

0.100 

0.100 

0.100 

0.100 

0.102 

0.110 

0.110 

0.110 

0.113 

0.118 

0.119 


-0. 
-0.037 
-0.100 
0.011 
0.039 
0.058 
0.087 
0.062 
0.085 
0.120 
0.142 
0.164 
0.180 
1.000 


1.000 

0.700 

0.250 

0.070 

-0.025 

-0.100 

-0.145 

-0.170 

-0.182 

-0.200 

-0.215 

-0.226 

-0.241 

-0.255 

-0.245 

-0.280 

-0.280 

-0.280 

-0.270 

-0.260 

-0.250 

-0.221 

-0.192 

-0.145 

-0.110 

-0.06l 

-0.030 

-0.005 

0.020 

0.050 

0.030 

0.065 

0.100 

0.142 

0.164 

0.180 

1.000 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


-0.110 

-0.114 

-0.117 

-0.119 

0.119 

0.116 

0.110 

0.100 

0.087 

0.070 

0.049 

0.025 

0.002 

0.029 

0.058 

0.086 

0.113 

0.137 

0.159 

0.176 

0.193 

1.000 


-1.000 

-0.178 

-0.166 

-0.167 

-0.170 

-0.173 

-0.178 

-0.184 

-0.190 

-0.198 

-0.207 

-0.217 

-0.228 

-0.239 

-0.249 

-0.258 

-0.256 

-0.269 

-0.270 

-0.266 

-0.258 

-0.244 

-0.225 

-0.200 

-0.169 

-0.135 

-0.097 

-0.057 

-0.017 

0.023 

0.06l 

0.096 

0.126 

0.152 

0.173 

0.191 

1.000 


TABLE  3 


Section  Shape  for  Duct  II 


1  -  z 

Design 

Calculated  from 

Calculated  from 

Experimental 

Theoretical  C 

l 

s(z) 

c(z) 

s(z) 

c(z) 

s(z) 

c(z) 

0 

0.00000 

0.00000 

0.0000 

0.0000 

0.0000 

0.0000 

0.05 

0.02095 

0.01085 

0.0185 

0.0075 

0.0210 

0.0125 

0.10 

0.02920 

0.01793 

0.0271 

0.0140 

0.0296 

0.0204 

0.15 

0.03528 

0.02348 

0.0338 

0.0193 

0.0359 

0.0254 

0.20 

0.0U002 

0.02797 

0.0393 

0.0240 

0.0408 

0.0289 

0.25 

0.01*364 

0.03162 

0.0435 

0.0281 

0.0446 

0.0315 

0.30 

0.04637 

0.03454 

0.0466 

0.0314 

0.0474 

0.0339 

0.35 

0.0U832 

0.03681 

0.0489 

0.0341 

0.0494 

0.0362 

0.1*0 

0.01*952 

0.03846 

0.0504 

0.0360 

0.0507 

0.0383 

0.1*5 

0.05000 

0.03952 

0.0513 

0.0373 

0.0512 

0.0400 

0.50 

0.01*962 

0.04000 

0.0514 

0.0379 

0.0510 

0.0413 

0.55 

0.04846 

0.03988 

0.0507 

0.0379 

0.0500 

0.0418 

o.6o 

0.04653 

0.03914 

0.0494 

0.0374 

0.0482 

0.0413 

0.65 

0.04383 

0.03774 

0.0474 

0.0359 

0.0456 

0.0396 

0.70 

0.04035 

0.03557 

0.0442 

0.0336 

0.0421 

0.0367 

0.75 

0.03612 

0.03248 

0.0402 

0.0303 

0.0377 

0.0325 

0.80 

0.03110 

0.02811 

0.0351 

0.0260 

0.0324 

0.0272 

0.85 

0.02532 

0.02170 

0.0289 

0.0207 

0.0262 

0.0209 

0.90 

0.01877 

0.01434 

0.0215 

0.0141 

0,0192 

o.oi4o 

0.95 

0.01143 

0.00  685 

0.0130 

0.0069 

0.0115 

0.0069 

1.00 

0.00000 

0.00000 

0.0000 

0.0000 

0.0000 

0.0000 

a  - 

0 

a  =  0.0048  degrees 

a  =  -0.0720  degrees 
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TABLE  4 


Section  Shape  for  the  BTZ  Duct 


0.00 

0.05 

0.10 

0.15 

0.20 

0.25 


0.30 

0.35 

0.40 

0.^5 

0.50 

0.55 

o.6o 

0.65 

0.70 

0.75 

0.80 

0.85 

0.90 

0.95 

1.00 


Design 

Calculated  from 

Calculated  from 

Experimental 

Theoretical  C 

P 

s(z) 

c(z) 

s(z) 

c(z) 

s(z) 

c(z) 

0.00000 

0.00000 

0.0000 

0.0000 

0.0000 

0.0000 

0.01257 

0.00000 

0.0114 

0.0002 

0.0129 

0.0001 

0.01752 

0.00000 

0.0172 

-0.0003 

0.0178 

0.0001 

0.02119 

0.00000 

0.0213 

-0.0006 

0.0214 

0.0002 

0.02401 

0.00000 

0.0245 

-0.0007 

0.0242 

0:0002 

0.02618 

0.00000 

0.0270 

-0.0007 

0.0265 

0.0002 

0.02782 

0.00000 

0.0288 

-0.0007 

0.0284 

0.0003 

0.02899 

0.00000 

0.0301 

-0.0008 

0.0298 

0.0003 

0.02971 

o.ooc^r 

0.0309 

-0.0008 

0.0307 

0.0004 

0.03000 

O.OOOoO 

0.0310 

-0.0007 

0.0311 

0.0004 

0.02985 

0.00000 

0.0304 

-0.0005 

0.0310 

0.C003 

0.02925 

0.00000 

0.0292 

-0.0003 

0.0302 

0.0003 

0.02815 

0.00000 

0.0273 

0.0000 

0.0287 

0.0C03 

0.02611 

0.00000 

0.0246 

0.0003 

0.0266 

0.0003 

0.02316 

0.00000 

0.0211 

0.0006 

0.0237 

0.0002 

0.01953 

0.00000 

0.0174 

0.0009 

0.0202 

0.0002 

0.01543 

0.00000 

0.0134 

0.0011 

0.0161 

0.0002 

0.01107 

0.00000 

0.0090 

0.0005 

0.0116 

0.0001 

0.00665 

O.OOOOO 

0.0052 

0.0002 

0.0070 

0.0001 

0.00262 

0.00000 

0.0023 

0.0000 

0.0027 

0.0000 

0.00000 

0.00000 

0.0000 

0.0001 

0,0000 

0.0000 

a  =  0 


a  =  0.0973  degrees 


a  =  -0.0097  degrees 


-.-a 


•K~ 


•X 

M 


■1 


% 


l 


'1 


percent  larger  than  the  design  thickness.  It  is  not  possible  to  draw  any 
general  conclusions  about  the  accuracy  of  predicting  the  thickness  except 
it  is  reasonable  and  probably  within  the  accuracy  of  either  the  linearized 
or  experimental  pressure  distributions.  This  is  all  that  can  be  expected. 

The  camber  appears  to  differ  by  a  somewhat  greater  magnitude  than 
the  thickness  for  Duct  II.  In  fact,  for  the  linearized  pressure  distri¬ 
bution  the  calculated  camber  is  somewhat  greater  than  the  actual  camber 
and  for  the  experimental  pressure  distribution,  the  camber  is  somewhat 
less.  The  section  angle  of  attack  is  less  than  the  design  angle  of  zero 
degrees  for  the  linearized  pressure  distribution  and  slightly  greater  for 
the  experimental  pressure  distribution.  Both  the  camber  and  angle  of 
attack  show  the  effect  of  accuracy  of  determining  the  pressure  distri¬ 
bution  near  the  leading  edge.  For  the  linearized  pressure  distribution, 
the  pressure  of  the  duct  is  calculated  to  be  infinite  at  the  leading  edge, 
whereas  for  the. experimental  pressure  distribution  the  closest  point 
measured  to  the  leading  edge  was  0.025  cf  the  chord.  In  either  case  the 
true  pressure  distribution  near  the  leading  edge  is  not  known. 

For  the  BTZ  duct  the  camber  and  angle  of  attack  calculated  from  the 
linearized  pressure  distribution  are  certainly  within  the  accuracy  that 
can  be  expected.  The  angle  of  attack  is  within  0.01  of  a  degree  and  the 
camber  is  only  1  percent  of  the  thickness  and  should  be  considered 
negligible.  Surprisingly,  the  calculated  camber  and  angle  of  attack  do 
not  appear  to  be  as  accurate  from  the  experimental  pressure  distribution 
as  from  the  linearized.  There  are  only  two  experimental  points  (one  in¬ 
side  and  the  other  outside)  within  10  percent  of  the  leading  edge  and  both 
the  calculated  results  of  Reference  7  and  the  present  calculations  indi¬ 
cate  that  the  point  on  the  inside  of  the  duct  is  incorrect. 

The  foregoing  results  of  the  ducts  by  themselves  are  presented  to 
give  an  indication  of  the  accuracy  of  the  procedure.  In  reality,  it  makes 
little  sense  to  use  a  linearized  theory  to  calculate  a  shape  from  a 
linearized  pressure  distribution.  The  usefulness  of  the  program  presented 
in  this  report  is  to  design  a  duct  for  a  given  pressure  distribution  when 
the  duct  is  in  the  presence  of  a  propeller  or  some  other  axisymmetric 
body.  This  allows  a  duct  shape  to  be  chosen  which  will  not  cavitate  or 
separate. 
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TABLE  5 


Pressure  Distribution  for  a  NACA  66-010  Thickness  Distribution 
and  Propeller  Induced  Velocities 


Leading  Edge 
0.0050 
0.0075 
0.0125 
0.0250 
'0.0500 
0.0750 
0.1000 
0.1500 
0.2000 
0.2500 
0.3000 
0.3500 
0.4000 
0-.  U500 
0.5000 
0.5500 
0.6000 
0.6500 
0.7000 
0.7500 
0.8000 
0.8500 
0.9000 
0.9500 
1.0000 


V 

c 

p~ 

*r>p/cT 

1.000 

1.000 

-0.0382 

0.104 

o.io4 

-0.0384 

0.028 

0.028 

-0.0386 

-0.023 

-0.023 

-0.0388 

-0.078 

-0.078 

-0.0389 

-0.125 

-0.125 

-0.0410 

-0.154 

-0.154 

-0.0426 

-0.174 

-0.174 

-0.0449 

-0.198 

-0.198 

-0.0480 

-0.215 

-0.215 

-0.0525 

-G.226 

-0.226 

-0.0564 

-0.236 

-0.236 

-0.0555 

-0.243 

-0.243 

-o.o4oo 

-0.249 

-0.249 

-0.0330 

-0.255 

-0.255 

-0.0210 

-0.261 

-0.261 

-0.0000 

-0.265 

-0.265 

0.0210 

-0.270 

-0.270 

0.0330 

-0.250 

-0.250 

o.o4oo 

-0.190 

-0.190 

0.0555 

-0.121 

-0.121 

0.0564 

-0.052 

-0.052 

0.0525 

0.021 

0.021 

o.o48o 

0.096 

0.096 

0.0449 

0.179 

0.179 

o.o4io 

0.271 

0.271 

0.0382 

t  _JL\  v,  /n 
V  T 


-0.0157 

-0.0158 

-0.0l60 

-0.0168 

-0.0176 

-0.0180 

-0.0199 

-0.0200 

-0.0230 

-0.0280 

-0.0697 

-0.1470 

-0.1660 

-0.1760 

-0.1860 

-0.1950 

-0.1860 

-0.1760 

-0.1660 

-0.1470 

-0.0697 

-0.0280 

-0.0230 

-0.0200 

-0.0180 

-0.0157 


w/o  Prop, 
w  Prop.  CT*I.O 
w  Prop.  Cr*2.0 


0.5 

1  -  Z 


Figure  5  -  Shape  of  Duct  with  Pressure  Distribution 
Corresponding  to  a  NACA  66-010  Thickness 
Form  with  and  without  a  Propeller 
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4.2  DUCTED  PROPELLER 


The  shape  of  the  duct  of  a  ducted  propeller  can  be  calculated  by 
the  program  by  inputting  the  steady  axial  and  radial  velocities  induced 
by  the  propeller  on  t2ie  duct.  These  velocities  for  use  in  this  manner  are 
discussed  in  References  9  and  14  and  have  been  tabulated  in  Reference  19. 
Calculations  for  duct  shape  were  made  for  a  pressure  distribution 

20 

corresponding  to  the  distribution  on  a  NACA  66-010  basic  thickness  form. 
Table  5.  Two  values  of  the  propeller  thrust  loading  coefficient  Cp  were 
assumed,  i.,e.,  Gp  =  1.0  and  Cp  =  2.0. 


where 

R  =  propeller  radius, 

T  =  propeller  thrust, 

V  =  velocity,  and 
p  =  mass  density. 

Propeller-induced  velocities  which  were  input  into  the  program  are  shown 
in  Table  5.  For  these  calculations  the  tip  clearance  was  assumed  to  be 
one  percent  of  the  propeller  radius  and  the  propeller  location  was 
assumed  to  be  at  the  midchord  of  the  duct.  Also,  for  convenience,  the 
thrust  loading  coefficient  was  based  on  the  average  velocity  at  the  pro¬ 
peller  and  not  the  free-stream  velocity  as  would  be  the  usual  case. 

Results  of  these  calculations  are  tabulated  in  Table  6  and  plotted 
in  Figure  5.  Also  plotted  in  this  figure  is  the  shape  of  the  duct  without 
the  propeller.  Not  only  is  the  camber  and  angle  of  attack  of  each  duct 
changed  considerably  but  the  thickness  is  also  changed.  The  biggest  effect 
is  on  the  angle  of  attack  which  goes  from  essentially  zero  fo  5.35  degrees 
for  a  Gp  =  1.0  and  from  zero  to  10.48  degrees  for  a  Cp  =  2.0.  Also  the 
cambers  of  all  the  ducts  are  S-shaped. 

With  the  propeller  in  the  duct,  the  duct  sections  are  thinner  near 
the  leading  edge  and  thicker  toward  the  trailing  edge  than  for  the  duct 
alone.  This  effect  increases  with  propeller  loading.  Also,  the  maximum  duct 
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TABLE  6 


Shape  of  Duct  with  and  without  a  Propeller 


0.00 

0.05 

0.10 

0.15 

0.20 

0.25 

0.30 

0.35 

0.1*0 

0.U5 

0.50 

0.55 

0.60 

0.65 


Without 

Propeller 

With  Propeller 

cT  =  1.0 

With  Propeller 

CT  =  2.0 

s(z) 

e(z) 

s(z) 

c(z) 

s(z) 

c(z) 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0192 

0.0020 

0.0133 

0.0054 

0.0073 

0.0088 

0.0281* 

0.0039 

0.0203 

0.0106 

0.0122 

0.0174 

0.0352 

0.0056 

0.0257 

0.0157 

0.0163 

0.0259 

o.oi*oi* 

0.0070 

0-.0302 

0.0207 

0.0201 

0.0343 

0.01*1*5 

0.0083 

0.0343 

0.0243 

0.0240 

0.0402 

0.01*77 

0.0093 

0.0382 

0.0244 

0.0287 

0.0394 

0.0500 

0.0101 

0,0422 

0.0222 

0.0344 

0.0343 

0.0511* 

0.0106 

0.0457 

0.0191 

0.0401 

0.0276 

0.0521 

0.0109 

0.0490 

0.0154 

0.0460 

0.0199 

0.0518 

0.0109 

0.0519 

0.0109 

0.0520 

0.0108 

0.0507 

0.0106 

0.0539 

0.0061 

0.0571 

0.0016 

0.0485 

0.0101 

0.0542 

0.0015 

0.0600 

-0.0070 

0.0445 

0.0092 

0.0525 

—0.0029 

o.o6o4 

-0.0151 

0.0389 

0.0081 

0.0485 

-0.0072 

0.0582 

-0.0225 

0.0322 

0.0068 

0.0426 

-0.0089 

0.0529 

-0.0247 

0.0250 

0.0054 

0.0353 

-0.0080 

0.0455 

-0.0214 

0.0177 

0.0039 

0.0272 

-0.0062 

0.0368 

-0.0164 

0.0106 

0.0025 

0.0188 

-0.0044 

0.0270 

-0.0112 

0.0044 

0.0012 

o.oio4 

-0.0022 

0.0164' 

-0.0056 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

a  =  0.13  degrees 


a  =  5-35  degrees 


a  =  10.1*8  degrees 


thickness  increases  with  propeller  loading.  Both  effects  are  a  consequence 
of  assuming  the  propeller  location  to  be  at  midchord. 

If  the  propeller  induces  a  velocity  in  the  same  direction  as  the 
free-stream  velocity  everywhere  on  the  duct,  the  thickness  is  decreased 
everywhere.  This  would  be  the  case  when  the  propeller- is  located  at  or 
aft  of  the  duct  trailing  edge.  Conversely,  if  the  propeller  induces  a 
velocity  in  the  opposite  direction  to  the  free-stream  velocity  everywhere 
on  the  duct,  the  thickness  is  increased  everywhere.  This  would  be  the 
case  when  the  propeller  is  located  at  or  forward  of  the  duct  leading  edge. 

The  calculations  for  a-  duct  shape  with  a  propeller  were  based  on 
the  assumption  that  the  inside  and  outside  pressures  were  the  same. 
Additional  calculations  were  made  for  a  thrust  loading  coefficient  of  two, 
Gj,  =  2.0,  and  the  pressure  distribution  given  in  Table  7.  The  pressure 
distributions  correspond  to  a  NACA  66-010  thickness  distribution  with  a 
NACA  a  =  0.8  mean  line  of  4-percent  camber.  In  one  case  the  camber  is 
negr.tJ>e  (toward  the  inside  of  the  duct)  and  in  the  other,  positive  which 
has  the  effect  of  shifting  the  negative  pressure  from  the  inside  to  the 
outside  of  the  duct. 

Results  of  these  calculations  are  shown  in  Table  8.  The  pressure 
distribution  with  a  negative  pressure  on  the  inside  of  the  duct  results  in 
a  decrease  in  duct  thickness  whereas  that  with  a  positive  pressure  on  the 
inside  of  the  duct  results  in  an  increase  in  duct  thickness.  This  is  a 

rWa(z)1 

consequence  of  the  change  in  sign  of  ^ — ^“^Jdiff  in  Ecluation  [2.2.16].  In 
fact,  the  effect  of  this  term  is  so  large  that  the  calculations  show  that 
to  achieve  the  pressure  distribution  chosen  for  the  calculation  with  a 
negative  pressure  inside  the  duct,  the  thickness  would  have  to  be  negative 
near  the  leading  edge.  Since  ducts  of  this  type  present  construction 
problems,  it  can  only  be  concluded  that  such  a  pressure  distribution  cannot 
be  achieved  for  the  propeller  loading  and  location  assumed. 

As  in  Table  6,  both  cambers  given  in  Table  8  are  somewhat  S-shaped. 
The  difference  being  that  for  the  negative  pressure  inside  the  duct,  the 
camber  is  generally  negative  and  for  the  positive  pressure  inside  the  duct, 
the  camber  is  generally  positive.  It  should  be  noted  that  the  angle  of 
attack  varies  less  than  a  degree  between  the  two  ducts. 
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TABLE  7 


Pressure  Distribution  of  a  NACA  60-OIO  Thickness  Form 
with  a  NACA  a  =  0.8  Mean  Line  of  h  Percent 


Positive  and  Negative  Camber 


]  -  2 

Positive  Camber 

negative 

Camber 

V 

C 

P“ 

V 

c 

p- 

Leading  Edge 

1.000 

1.000 

1.000 

1.000 

0.0050 

-0.23*1 

0.387 

0.387 

-0.231* 

0.0075 

-0.323 

0.321* 

0.321* 

-0.323 

0.0125 

-0.378 

0.283 

0.283 

-0.378 

0.0250 

-0. St  1(5 

0.236 

0.236 

-0. 1*1*5 

0.0500 

-0.501 

0.195 

0.195 

-0.501 

0.0750 

-0.533 

0.172 

0.172 

-0.533 

0.1000 

-0.558 

0.151* 

0.151* 

-0.558 

0.1500 

-0.585 

0.133 

0.133 

-0.585 

0.2000 

-0.603 

0.120 

0.120 

-0.603 

0.2500 

-0.615 

0.111 

0.111 

-0.615 

0.3000 

-0.628 

0.101 

0.101 

-0.628 

0.3500 

-0.636 

0.096 

0.096 

-0.636 

0.1*000 

-0.61*1* 

0.090 

0.090 

-0.61*1* 

0.1*500 

-0.61*9 

0.086 

0.086 

-0.61*9 

0.5000 

-0.656 

0.080 

0.080 

-0.656 

0.5500 

-0.662 

0.076 

0.076 

-0.662 

0.6000 

-0.667 

0.073 

0.073 

-0.667 

0.6500 

-0.61*1* 

0.090 

0.090 

-0.61*1* 

0.7000 

-0.575 

0.11*1 

0.11*1 

-0.575 

0.7500 

-0.1*96 

0.199 

0.199 

-0.1*96 

0.8000 

-0.1*16 

0.257 

0.257 

-0.1*16 

0.8500 

-0.237 

0.250 

0.250 

-0.237 

0.9000 

-0.067 

0.21*5 

0.21*5 

-0.067 

0.9500 

-0.103 

0.252 

0.252 

-0.103 

1.0000 

1.000 

1.000 

1.000 

1.000 

TABLE  8 


Shape  of  Duct  with  a  Propeller  for  Positive 
and  Negative  Inside  Pressures 


1  -  z 

From  Fressure  Distribution 
with  Positive  Camber 

From  Pressure  Distribution 
with  negative  Camber 

s(z) 

c<z) 

s(z) 

c(z) 

0.00 

0.0000 

0.0000 

0.0000 

0.0000 

0.05 

0.0292 

0.0212 

-0.0086 

-0.0030 

0.10 

0.01*28 

0.0381 

-0.0103 

-0.0022 

0.15 

0.0532 

0.0531 

-0.0110 

0.0001 

0.20 

0.0619 

0.0668 

-0.0110 

0.0037 

0.25 

0.0697 

0.0770 

-0.0100 

0.0056 

0.30 

0.0771* 

0.0795 

-0.0075 

0.0016 

0.35 

0.0853 

0.0771 

-0.0035 

-0.0059 

0.1*0 

0.0921* 

0.0723 

0.0010 

-0.011*3 

0.1*5 

0.C991 

0.0656 

0.0063 

-0.0231 

0.50 

6.1051* 

0.0569 

0.0122 

-0.0325 

0.55 

0.1100 

0.01*71* 

0.0177 

-C.0l*l>< 

0.60 

0.1118 

0,0376 

0.0215 

-0.01*89 

0.65 

O.llOl* 

0.0275 

0.0233 

-0.0551 

0.70 

0.1056 

0.0170 

0.0230 

-0.0597 

0.75 

0.0972 

0.0107 

0.0203 

-0.0580 

0.80 

0.0857 

0.0085 

0.0l6l 

-0.01*96 

0.85 

0.0717 

0.0062 

0.0110 

-0.0376 

0.90 

O.C561 

0.0037 

0.0061* 

-0.0251 

0.95 

0.0378 

0.0027 

0.0027 

-0.0132 

1.00 

0.0000 

0.0000 

0.0000 

0.0000 

a  =  11.01  degrees 

a  =  9-95  degrees 
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CONCLUSIONS 


This  report  presents  a  computer  program  for  the  inverse  problem  of 
the  annular  airfoil.  As  a  result  of  calculations  made  with  this  program 
the  following  conclusions  can  be  made: 

1.  With  the  restrictions  of  linearized  theory,  duct  shapes  can  be 
determined  for  desired  pressure  distributions  even  in  the  presence  of  a 
propeller. 

2.  For  a  given  pressure  distribution,  the  presence  of  the  propeller 
at  the  duct  midchord  tends  to  increase  the  section  angle  of  attack, 

make  the  camber  S-shaped,  and  move  the  point  of  maximum  thickness  toward 
the  trailing  edge. 

3.  For  a  positive  pressure  on  the  inside  of  the  duct,  the  duct 
thickness  is  increased.  For  a  negative  pressure  on  the  inside  of  the  duct 
and  a  positive  pressure  outside  the  duct  thickness  is  decreased. 

4.  The  propeller  location  and  loading  have  important  effects  on  the 
duct  shape. 

5.  The  computer  program  is  quite  versatile  and  will  facilitate  the 
design  of  ducted  systems. 
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APPENDIX  A  -  Continued 
OUTPUT 


THE  FOLLOWING  TABLE  IS  THE  DISTRIBUTION  OF  THE  DESIRED  FOIL 
AT  A  5.34884  DEGREE  ANGLE  OF  ATTACK. 

ABSCISSA  THICKNESS  CAMBER 


leading  edge 
o.os 
0.10 
C.1S 
0.20 
0.25 
0.30 
0.35 
0.40 
0.45 
0.50 
0.55 
0.60 
C  .65 
0.70 
0  .75 
0.80 
0*85 
0.90 
0.95 

TRAILING  EDGE 


0. 

C. C 1 3274 
0.020276 
0.C25718 
0.030249 
C. 034286 
0.038213 
0.042190 
0.045733 
0.049012 
0.051928 
0.CS3906 
0.054231 
0.052961 
0.048536 
0.042579 
0.035285 
0.027224 
0.018823 
0.010426 
0.000053 


0. 

0.005393 

0.010640 

0.015720 

0.020667 

0.024266 

0.024375 

0.022209 

0.019137 

0.015378 

0.010850 

0.006127 

0.001526 

-0.002947 

-0.007205 

-0.008922 

-0.008012 

-0.006227 

-0.00435S 

-0.002209 

-0.000003 


APPENDIX  B  -  FORTRAN  LISTING  OF  COMPUTER  PROGRAM 


.  t 
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£ 


o 
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FOR  YV18 

DIMENSION  POUT  !37)'»PIN(37)  .PAX ( 37 ) .PRAD! 37 ) . THETA! 37 > .THETAB!  101 ) 

1  *  THICK  ( 101 )  .OMEGA ( 37  )  .CPRIME  ( 101  >  ♦C(21)»S(21)»Z(21)*X(3N7) 

2  »DUM (21 ) 

READ  199  jNCASE 
199  FORMAT  (14) 

DO  198  ICASE  *  l.NCASE 
READ  99 

99  FORMAT  (72K1 
1 

READ  10,  H.WX0.M.N.NZ.NP.NH.IJK.J1.IPUN 

10  FORMAT  ( 2F10.6  #8 14 ) 

READ  11#  (X(I>.POUT(I).PlN(I).PAX(I).PRAD(n.  I=1»NZ) 

11  FORMAT  (5F14.8) 


WAKE  -  1.0  -  WXD  YV18 

PRINT  99 

PRINT  13#  H.WAKE  YV18 

13  FORMAT  (98H0  CALCULATION  OF  THE  THICKNESS  DISTRIBUTION  AND  THE  CAM 
1BER  DISTRIBUTION  USING  THE  FOLLOWING  DATA.  //22H  CHORD  DIAMETER  R 
2AT-I0  5X.15H  WAKE  FRACTION  /7X#F12.8  *9X#F12.8//17H  ABSCISSA  OF  TK 
3E  10X.16H  OUTER  PRESSURE  5X.16H  INNER  PRESSURE  10X.24H  INDUCED  VEL 
40CITY  TERMS  /23H  PRESSURE  DISTRIBUTION  52X.7H  AXIAL  6X»8H  RADIAL 


5  //) 

PRINT  14#  (X  <I) »POUT(I ) .PIN! 
14  FORMAT  (1H  8X  .F8.5 .1-4X.F12.7 .9X 
DO  765  I  =  1  .N 
765  PRAD(I)  =>  -PRAD! I) 

PI  =  3.14159265 
ON  *  N-l 
DELTA  *  PI/DN 
NJN  =  N-3 
DO  115  I = 1 #NZ 
X(I)  ■  1,0  -  X ( I ) 

XARG  =  2.0*X( I )  -  1.0 
THETA! I )  =  ACOSF ! ABSF ( XARG )  ) 

IF  (XARG)  114. 115.115 

114  THETA! I )  =  PI  -  THETA! I) 

115  CONTINUE 
THETAB(l)  =  0.0 


) .PAX! I ) .PRAD! 1 ) . 1=1 *NZ ) 
F12.7.12X.F12.7.F14.7) 


YV18 

YV18 

YV18 

YV13 


DO  116  1=2. N 

116  THETAB(I)  =  THETAP(I-l)  +  DELTA 
18  IF  ( ABSF ( POUT ( 1 )  -PIN(l))  -  .0001)  210.210.211 

210  fZ  *  1 

GO  TO  212 

211  IZ  *  2 


10 


80 

90 


180 

190 

200 

220 


212  DO  25  I=IZ.NZ 

OMEGA! I )  =  PI*(.25*(P0UT< I)  +  P IN ! I > )  -  PAX  (I)) 

25  PAX  (I)  =  . 5* { POUT ( I )  -  PIN! I  I >  *S INF ( THETA! I ) ) 

C  PAX  IS  NOW  THE  VELOCITY  DIFFERENCE  TERM  TIMES  SINE!  THETA  ). 
GO  TO  (27  .225  )  .  IZ 
225  OMEGA! 1 )  =  -PI»PAX(1> 

DO  224  1  =  I Z »NZ 
THETA! 1-1 )  =  THETA! I ) 

224  PAX(I-l)  =  PAX  !  I ) 

NZ1  =  NZ-1 

CALL  DISCOT  (0.0.0, 0.THETA.PAX.PAX.-020.NZ1.0.PAXT) 

DO  231  1=1 »NZ 
POUT! I )  =  THETA! I ) 

231  PIN! I )  =  PAX ( I ) 

DO  230  I  =2  »NZ 


36 


THETA(I)  -  POUT  (1-1 ) 

230  PAX ( I )  =  P IN ( 1-1 ) 

THETA(l)  *  0„0 
PAX ( 1 )  =  PAXT 
27  Z<1>  =  0.0 
DO  30  I  * 2 » 2 1 
30  Z (I)  *  Z(I-l)  +  .05 

CALL  THKDIS  (DELTA. DN.  H.M.N.NJN.NZ.PI .THICK.  YV18  230 

1 PAX. OMEGA .NP.NH  . THETA .THETAB  »Z.S» I JK> 

CALL  CAMBER  (DELTA. DN.  H.  N .NJN.NZ *P I »TH ICK .WAKE.PAX.PRAD. J1 » THETA 
IB f THETA. ATTACK .CPR I  ME. S.Z) 

CALL  SOLVE  ( I JK.P ( .THETAB.CPRIME .N.Z .C ) 

PRINT  98.  ATTACK 

98  FORMAT  (67H1  THE  FOLLOWING  TABLE  IS  THE  DISTRIBUTION  OF  THE  DESIRE 
ID  FOIL  AT  A  F9.5.25H  DEGREE  ANGLE  OF  ATTACK.  /15X.9H  ABSCISSA  1 
20X.10H  THICKNESS  8X.8H  CAMBER  //) 

PRINT  150.  S(21 ) »C(21 ) 

150  FORMAT  (11X.14H  LEADING  EDGE  9X.F10.6  .F16.6 ) 

DO  777  1-2.20 

JIP  =  22*1 

777  PRINT  97.  Z  CD  .S  ( JI P )  »C(  JI P  ) 

97  FORMAT  ( 10X . F12 . 2 , 12X . F10. 6 »F 16.6 ) 

PRINT  151»S(1 J»C(1) 

151  FORMAT ( 10X. 15H  TRAILING  EDGE  9X .F] 0.6 »F 16.6 J 
IF  ( I PUN  -  II  198.198.599 

599  PUNCH  600.  (Z (  I )  »  1  =  1 .21 ) 

DO  610  1=1.21 

JIP  =  22-  I 

610  DUM(I)  =  C(JIPJ 

PUNCH  600.  (DUM(I). 1=1.21) 

DO  611  1=1.21 
JIP  =  22-  I 

611  DUM(I)  =  S(JIP) 

PUNCH  600.  (Z(  11.1  =  1.21) 

PUNCH  600.  ( DUM (I). 1=1*21) 

600  FORMAT  (9F8.5) 

198  CONTINUE 

CALL  ENDJOB  YV18  310 

END  YV18  320 

FOR  YV181010 

SUBROUTINE  THKDIS  ( DELTA .DN .H »M »N »NJN .NZ »P I » TH I CK »PAX .OMEGA »NP »NH 
1 .THETA . THETAB. Z.S. I JK) 

DIMENSION  ALPHA! 4C» 40) .B (40.101 ) .GAMMA ( 40 .40 ) . FM( 40 ) . FTHETB ( 10 1 ) . 

10MEGA( 37) .THETAB (103 ) .THICK! 101 ) .PAX ( 37 ) ♦ THETA ( 37 )  . 


EQUIVALENCE 

ERROR  =  0.00001 

CALL  KERNEL  ( B. DELTA. H.M . 


(C'»MMA  .ALPHA) 


CALL  KERNEL  ( B. DELTA. H.M .  N. PI .ERROR .KEEP ) 

LOOP  FOR  CALCULATING  A  SERIES  OF  F( THETA  BAR)  VALUES  BY  INCREMENTING 
THETA  BAR  FROM  0.0  TO  P I . 

DO  5?  1=1 »N 
IF  (1-1)  52.52.51 

51  IF ( I -N ) 54.52 ♦ 52 

54  SINB  =  1.0/(2.0*P.T#SINF(  THETAB (  I  )  )  ) 

52  CALL  USEGM  (THETAB! I ) »SINB» I .PI .NP.N.DN.H. THE TA.PAX.NZ .OMEGA >NH. 
1FTHETB ( I ) ) 

CALL  CALCAL  ( ALPHA, B.DELTA.  KEEP ,N . THETAB. P I 

CALL  CALCFM  ( B. DELTA. FM » FTHETB .  KEEP.N  ) 

LOOPS  FOR  PREPARATION  OF  COEFFICIENTS  OF  LINEAR  SYSTEM  FOR  MAT  I NV 
61  DO  70  I =1 .KEEP 
DO  70  J=1 .KEEP 


YV070250 

YV070260 

YV070270 

YV070290 

YV070315 

YV070330 


YV070430 

YV070440 

YV070520 

YV070530 

YV070530 
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-T* 


1 

,i 


$  .j 


i 

5 


! 


! 

f 

j 


I 


IF  (I-J)  71*72*71 

72  GAMMA(IjJ)  =  1.0  -  ALPHA (1  * J  5 /P  I 
GO  TO  70 

71  GAMMA! I  * J  )  =  -ALPHA ( I ♦ J ) /P I 
70  CONTINUE 

CALL  MAT  1 NV  ( GAMMA. KEEP » FM ♦ 1 *X *  I D ) 

I F { 2-ID)  75*76*75 

76  PRINT  77 

77  FORMAT  ( 1 17H0COEFF IC I  ENT  MATRIX  IS  SINGULAR. 


YV070540 

YV070560 

YV070570 

YV070590 

THE  INTEGRAL  EQUAT I0YV070600 


IN  IS  EITHER  INSOLVABLE  OR  HAS  AN  INFINITE  NUMBER  OF  SOLUTIONS. 
CALL  ENDJOB 

C  LOOPS  FOR  CALCULATION  OF  THE  SLOPE  OF  THE  THICKNESS  DISTRIBUTION. 
75  DO  90  I = 1 *N 

IF  (1-1)  92*92.85 
85  IF(I-N)  86*93*93 

92  THICK  U)  =  FTHETC(l)  +  FM(1)/PI 
DO  150  J  =  2 *KEEP 

DJ  =  J 

150  THICK  (I)  *  THICK  (I)  +  FM(J)*DJ/PI 
GO  TO  90 

93  THICK  (I)  =  FTHETBfl)  +  FM(1)/PI 
DO  149  J  =  2. KEEP 

DJ  =  J 

149  THICK  (I)  =  THICK  (1)  -  (-1.0)**  J*FM ( J ) *D J/P 1 


GO  TO  90 

86  SINT=  SINF(THETABM  5  ) 

SUM  =  FM( 1 ) 

DO  91  J=2»KEEP 
DJ  =  J 

91  SUM  =  SUM  +  FM(J)*SINF(DJ*THETAB( I ) ) /SINT 
THICK  (I)  =  FTHETB(I)  +  SUM/PI 
90  CONTINUE 

CALL  SOLVE  ( I JK ,P I , THET AB, TH I CK  *N,2*S) 
RETURN 
END 
FOR 


SUBROUTINE  KERNEL  ( B.DELTA*H»M*N  *P I  *  ERROR » KEEP ) 

DIMENSION  B( 40 • 101 ) »CONSK{ 101 ) *EE ( 10 1 ) . FB ( 101 ) .THETA1I 101 ) * 
1THET  A5 { 201 ) *FBCOS(201) 


C  LOOP  FOR  VARYING  THETA  1  IN  CALCULATION  OF  FOURIER  COEFFICIENTS  FOR 
C  APPROXIMATION  OF  THE  KERNEL. 

KEEP  =  1 
DO  200  I  *  1 *N 
IN  =  1 

C  LOOP  FOR  CALCULATION  OF  ACTUAL  FOURIER  COEFFICIENTS  —  8(THETA  1). 


YVO7O6IO 

YV070700 


YV070720 

YV070730 

YV070740 

YV070750 

YV070760 


YV181080 
YV181090 
YV14  10 

YV140030 

YV14C040 

YV140050 

YV140060 

YV140070 

YV140080 

YV140090 


DO  202  L=1*M 

YV140100 

DL  =  L 

YV140110 

11  =  1 

YV140120 

IF  (1-1)  700*700*710 

YV140130 

700  THETAl(l)  =  0.0 

YV140140 

THETA1  VARIES  FROM  O.C  TO  PI* 

YV140150 

GO  TO  755 

YV140160 

710  THETAl(I)  =  THETAKI-l)  +  DELTA 

YV140170 

755  COS1  =  COSF ( THETA1 ( I ) ) 

YV140180 

IF  (10-L!  740*741,741 

YV140190 

741  GO  TO  (745*760)*  IN 

YV140200 

745  IN  =  2 

YV140210 

NN  *  51 

YV140220 

NNN  =  51 

YV140230 

ANGLE  =  PI/50.0 

YV140240 
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FIRST  =  0.0 
JJJ  =  1 
GO  TO  750 

740  IF  120-L)  742.743.743 
743  GO  TO  (746.746.760) .IN 

746  IN  =  3 

DO  201  J=1 .50 
LESS  =  103  -  2*J 
MINUS  =  52-J 

THETA5 ( LESS )  =  THETA5 (MI NUS ) 

201  FB(LESS)  =  FB ( MINUS ) 

NN  =  100 

NNN  =  101 
ANGLE  =  PI/100.0 
FIRST  =  ANGLE 
JJJ  =  2 
GO  TO  750 

742  GO  TO  (747.747.747.760) . IN 

747  IN  *  4 

DO  205  J=i»l00 

LESS  =  203  -  2*J 

MINUS  =  102  -  J 

THET  A5 ( LESS )  =  THLTA5(MINUS) 

205  FB (LESS )  =  FB(MINUS) 

NN  *  200 
NNN  =  201 
ANGLE  =  PI/200.0 
FIRST  =  ANGLE 
JJJ  »  2 

750  CALL  INTGRL  (ANGLE  .CONSK.COS1 .EE.FB. FIRST. I . 1 1 »L»  NN.NNN. 

1THET.M  (  I  )  .THETA5.JJJ.M.H) 

760  DO  203  J* 1 .NNN 

203  FBCOS(J)  =  F8( J)*COSF(DL*THETA5( J)) 

C  INTEGRAND  VALUES  FOR  INTEGRATION  WITH  RESPECT  TO  THETA5 
SUMS  =  FBCOS(l)  +  4. 0*FBCOS ( NNN-1 )  +  FBCOS(NNN) 

NJN  =  NNN-3 
DO  204  J=2 .NJN .2 

204  SUMB  =  SUMS  +  4.0*FBCOS(  J)  +  2 .0*FBC0S  ( J+l ) 

B(L.I)  =  2«0*ANGLE*SUM6/(3.0*PI  ) 

C  B  IS  THE  VALUE  OF  THE  FOURIER  COEFFICIENT  FOR  EACH  VALUE  OF  THETA  1. 

IF  ( ABSF ( B ( L . I ) )  -  ERROR)  790.790.202 
C  CHECK  FOR  THE  CONVERGENCE  OF  THE  FOURIER  COEFFICIENTS 

202  CONTINUE 
KEEP  =  M 
GO  TO  200 

790  L  =  L+l 

DO  771  LL=L.M 
771  B(LL.I)  =  0.0 
KEEP  1  =  L  -  1 

IF  (KEEP  -  KEEP1)  795.200.200 
795  KEEP  =  KEEP! 

200  CONTINUE 
RETURN 
END 


YV140250 
YV140260 
YV140280 
YV140290 
YV140300 
YV14031 0 
YV140320 
YV14  325 
YV14  326 
YV140330 
YV140340 
YV140350 
YV140360 
YV140370 
YV140380 
YV140400 
YV140410 
YV140420 
YV140430 
YV140440 
YV14  445 
YV14  446 
YV140450 
YV140460 
YV14C470 
YV14C480 
YV140490 
YV140500 
YV140520 
YV140530 

YV140560 

YV140570 

YV140580 

YV140590 

YV140600 

YV140610 

YV140620 

YV140670 

YV140680 

YV140690 

YV140700 

YV140710 

YV140770 

YV140840 

YV140850 


YV140860 

YV140870 

YV140880 


FOR 

SUBROUTINE  USEGM 


(THETAB.SINB.I .PI.NP.N.DN. H.THETA.WDIFFR. NZ. OMEGA 


l.NH.SUM) 

DIMENSION  THETA4(l0l).  WDIFFR(37  ) .OMEGA( 37  ! .ANSWER ( Id ) .  YV21  O3O 

1T( 101 ) »F(101 ) .R(507) »A( 100) .THETA! 37)  YV21  040 

IF  (1-1)  501.501.500  YV21  050 
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YV21  060 


501  LIFT  =  1 
NP1  =  NP/2 
NP2  *  NP1  +  1 
P  *  flP  1 
DELANG  =  PI/P 
THETA4 ( 1 )  =  0.0 
ANSWERIlt  =  0.0 

C  LOOP  FOR  CALCULATING  SERIES  OF  FUNCTIONAL  VALUES  FOR  HARMONIC  ANALYS 
00  510  J  =  2»NP2 

THETA4U)  =  THETAA(J-l)  +  DELANG 
TLOW  «  -CUBER  TF { P J  -  THETA4 ( J ) ) 

TUP  *  CUBER TF ( THFTA4( J) ) 

DELTAT  =  I  TUP  -  TLOW) /ON 
SIN2  =  SlNF ( THETA4 ( J ) )  #*  2 

C  LOOP  FOR  EVALUATION  OF  INTEGRATION  WITH  RESPECT  TO  T. 

DO  515  K=1*N 

IF  ( K-l )  516*516*517 

516  T ( 1 )  =  TLOW 
GO  TO  555 

517  TOC)  *  T t K-l  »  +  DELTAT 

555  IF  (ABSF(T(K)J  -  C.0001 )  556*556*55? 

556  FOC)  =  0.0 
GO  TO  515 

557  X  -  THETA4  ( J)  -  T(K)**3 
COST  =  COSF ( THETA4 ( J) ) 

COSX  =  COSF ( X ) 

C  CALCULATION  OF  ELLIPTIC  CONSTANT  K 

CONK  *  SORTFI l.O/I .25*(H*(COST  -  C0SX))**2  +  1.0)) 

C  SUBROUTINE  FOR  CALCULATION  OF  ELLIPTIC  INTEGRALS. 

CALL  WORK  (CONK. THETA4IJ). T(K).H,EK.EE. DIF) 

C  SUBROUTINE  FOR  INTERPOLATION  OF  VELOCITY  DIFFERENCE  VALUES. 

CALL  DISCOT  ( X *X . THET A ,WD I FFR *WD I FFR ♦ -1 20 *NZ .0 » WD I FR ) 

C  CALCULATION  OF  INTEGRAND  VALUES  FOR  INTEGRATION  WITH  RESPECT  TO  T. 

F  (  K  )  =  H*WDIFR*CONK*DIF»3«.0*T(K)»*2  /4.0 

515  CONTINUE 

SUMT  =  F(l)  +  4.0-F (N— 1 )  +  F(N) 

NTN  =  N-3 
DO  560  L=2*NTN»2 

560  SUMT  =  SUMT  +  4.0*F(L)  +  2.0*F(L+1) 

C  SUBROUTINE  FOR  INTERPOLATION  OF  MEAN  VELOCITY  MINUS  AVERAGE  INDUCED 
C  VELOCITY  VALUES 

CALL  DISCOT  ( THETA4 ( J )» THETA4 ( J )♦ THET A .OMEGA *OMEGA ♦ -120 »NZ *0* 
1Z0MECA) 

510  ANSWER ( J )  =  (DELTAT*SWMT/3.0  +  ZOMEGA)»S2N2 
DO  520  J=2.NP1 
NMINUS  =  NP  -  J  +  2 
520  ANSWER (NMINUS)  =  ANSWER (J) 

C  SUBROUTINE  FOR  CALCULATING  FOURIER  COEFFICIENTS. 

CALL  GMHAS  (NP ,NH, ANSWER *R ( 507 ) ) 

MAX  =  507 
DO  580  L=1*NH 
IL  =  MAX-5*L 
580  A(L)  =  R ( I L  ) 

500  SUM  =  A ( 1 ) / ( 2 • 0*P I ) 

LEFT  =  L I  FT 

GO  TO  *'570*571  )  »LEFT 
570  DO  572  L=2*NH 
OL  =  L 

572  SUM  =  SUM  +  DL*A ( L )/( 2 .0*P I ) 

LIFT  =  2 


YV21 
YV21 
YV21 
YV2 1 
ISYV21 
YV2  1 
YV2 1 
YV21 
YV21 
YV2 1 
YV2 1 
YV21 
YV2  1 
YV2 1 
YV21 
YV21 
YV21 
YV21 
YV2 1 
YV2 1 
YV2 1 
YV21 
YV21 
YV21 
YV21 
YV21 
YV21 
YV21 
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1 


GO  TO  599 

1 

571 

IF  ( I  — N )  575*574*574 

YV2 1  700 

575 

DO  576  L=2*NH 

YV21  710 

MM- 

DL  =  L 

YV2 1  720 

576 

SUM  =  SUM  +  A ( L ) *S INF ( DL*THETAB ) *S INB 

YV21  730 

GO  TO  599 

YV21  740 

' 

574 

DO  577  L=2  *NH 

YV2 1  750 

DL  =  L 

YV21  760 

577 

SUM  =  SUM  -  (-1.0)**L*DL*A(L) /(2.0*PI ) 

YV21  770 

599 

RETURN 

YV2 1  780 

END 

YV2 1  790 

FOR 

YV070010 

SUBROUTINE  CALCAL (ALPHA.B. DELTA,  KEEP *N * THSTA1 * P I 

YV070020 

DIMENSION  ALPHA(4C. 40) *0(40*101 ) *THETA1(101> 

YV070030 

NN=(N-3)/2 

YV070040 

DO  3:2  LL=i tKEbP 

DO  350  L=1 .KEEP 

YV070100 

DL  =  L 

YV070110 

SUMAL  =  B ( LL . 1 ) *DL  +  4 .0*8 ( LL *N- 1 ) *S I NF ( DL*THE T A1 ( N~ 1 ) ) /S I NF ( THETAYV070120 

ll(N-l))-  B(LL*N)*DL*COSF(DL*PI  ) 

DO  351  J=1 *NN  YV070150 


351  SUMAL  =  SUMAl.  +  4.0*B ( LL .2* J ) *S I NF ( DL*THETA1 ( 2* J ) ) /S 1 NF ( THETA1 ( 2* JYV070160 

D)  +  2.0*B(LL*2*J<  1)*SINF(DL*THETA1(2*J+1)  >/S!NF(THETAl(2*J+l>  )  YV070170 

350  ALPHA(LL.L)  =  DEL T A* SUMAL/ 3 o0  YV070180 

352  CONTINUE  YV070260 

RETURN  YV070270 

END  YV070280 

FOR  YV070010 

SUBROUTINE  CALCFM < B .DELTA, FM ,FTHET 1  *  KEEP.N  YV07 

DIMENSION  B(40.101 ) ,FM(40) .FTHETH 101 )  YV070030 

NN=(N-3)/2  YV070040 

DO  301  L=1.KEEP  YV070050 

SUMFM=B(L*1 )*FTHET1(1 )+4.0*6(L*N-l >*F  CHET  1 ( N-l ) +B ( L  *N ) *FTHET1 ( N)  YV070100 


DO  300  1=1, NN 

300  SUMFM=SUMFM+4.0*B(L»2*I > »FTHET] <2*1 ) +2 . 0*e ( L ♦ 2* I +1 > «FTHET 1 ( 2* I +1 )  YVO7O12O 

301  FM(L )  =  DELTA*SUMFM/3.0  YV070130 


RETURN 

END 

FOR 

SUBROUTINE  CHEAT  ( THETA4 t T ,H . D I F ) 

CON  =  H*(T**6/2.0*COSF(THETA4)  -  T**3*5INF(THETA4) ) 
CK  =  CON/SQRTF(4.0+CON**2) 

IF  (CK)  930*931 *930 

930  ALAM  =  LOGF ( ABSF(4 .O/CK )  ) 

DIF  =  ALAM  -  1.0  -  ALAM#CK# *2/4.0 
RETURN 

931  DIF  =  0.0 
RETURN 
END 

FOR 

SUBROUTINE  WORK  ( CONK  * THETA4 *T *H ,EK *  EE -  Dl F ) 

IF  (CONK  ~  0,995)  940*940*941 

940  CALL  ELLIP  ( 0 . 0 .CONK  * 7 « ZZ »EK » EE ) 

GO  TO  981 

941  IF  (CONK  -  0.999999)  *942.943 

942  CALL  VELLIP  ( CONK*EE»EK ) 
voi  \)  1  r  -  c<  -  cc 

RETURN 

943  CALL  CHEAT  ( THEf A4 *T *H*D IF ) 

EK  =  OIF  +  1.0 


YV070210 

YV07 

010 

YV07 

020 

YV07 

030 

YV07 

040 

YV07 

050 

YV07 

060 

YV07 

070 

YV07 

080 

YV07 

090 

YV07 

100 

YV07 

005 

YV07 

010 

YV07 

020 

YV07 

030 

YV07 

040 

YV07 

050 

YV07 

060 

VWA  1 

1  VU  * 

V  1  V 

YV07 

080 

YV07 

090 

YV07 

100 
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EE  =  -1.0 
945  RETURN 
END 
FOR 

SUBROUTINE  EELLIP  ICONSK.EE ) 

IF  (CONSK  -  0.995)  940.940.941 

940  CALL  ELLIP  i 0. 0. CONSK *Z .ZZ.ZZZ.EE) 

RETURN 

941  YK  =  1.0  -  C0NCK**2 
IF  { YK )  94?, 942.941 

941  ALAM  =  LOGF ( 4 .C/SQRTF ( YK )  ) 

EE  =  1.0+!2»0*ALAM-1.0)*YK/4.0  +  3.0* ( AL AM-1 3 .0/ 1 2 . 0 ) *YK**2/ 1 6.0 
RETURN 

942  EE  =  1.0 
RETURN 
END 

LOAD  BATCH 
FOR 

SUBROUTINE  VELLIP  (XK.EE. EK) 

C  EVALUATION  OF  COMPLETE  ELLIPTIC  INTEGRALS  WHEN  K  SQUARED  IS  GREATER 
C  THAN  .99  BUT  NOT  EQUAL  TO  1.0  BY  ASYMPTOTIC  SERIES. 

YK  *  1.0  -  XK**2 

ALAM  =  LOGF ( 4 .O/SQRTF ( YK ) ) 

EE  ■=  1.0  +  (2.0*ALAM-1.0)*YK/4.0  +  3  »  0* ( A LAM- 13.0/ 12.0 )*YK**2/16< 
EK  =  ALAM  +  < ALAM-1.0)*YK/4.0  +  9.0* ( ALAM-7 .0/6 .0 )*YK*« 2/64.0 
RETURN 
END 
FOR 

SUBROUTINE  SOLVE ( JK.PI ♦ THETAB »F ,N »Z ♦ STORE ) 

DIMENSION  7 HE TAB  1101) »F ( 101 ) .STORE! 2 1 ) .Z ( 2 1 ) 

I JK  =  JK 

DO  101  1  =  1  ,20 
DIJK=IJK-1 
ZARG=2.0*Z( I J-1.0 
ANG=ACOSF(ABSF(ZAPG) ) 

IF(ZARG)134. 115,135 

134  ANG=PI -ANG 

135  DEL=ANG/DIJK 

CALL  D I  SCOT ( DEL, DEL, THETAB.F.F  ,“120, N.O .RESULT) 

STORE! I ) =F ( 1 ) +4. 0*RESULT 

! JK= I JK-1 

DO  100  J=3  » I JK ,2 

D  J  =  J-l 

T=DJ*DEL 

CALL  DI  SCOT(  T  .T.THETAB » F  »F. --120  »N»0»  RESULT) 

STORE! I )=STORE! I )+2.0*RESULT 
T=(DJ+1.0i*DEL 

CALL  D I  SCOT! T ♦ T, THETAB. F, F .-1 20, N.Ot RESULT) 

100  STORE! I )=STORE( I )+4.0*RESULT 

CALL  D I  SCOT ( ANC, ANG, THET A8, F, F, - 120, N.O. RESULT) 

STORE! 1 ) = ( STORE { I )+RESULT)*DEL/3.0 

101  I JK= I JK-1 
STORE! 2 1 ) =0.0 
RETURN 

END 

FOR 

SUBROUTINE  CWORK  (CONK .TERM, EK »EE, DIF ) 

IF  (CONK  -  0,995)  940.940,941 
940  CALL  ELLIP  (0.0, CONK »ZZ >ZZZ ,CX »EE ) 

GO  TO  981 


YV07 
YV07 
YV07 
YV07 
YV07 
YVC7 
YV07 
YV07 
YV07 
YVO? 
YV07 
YV07 
YV07 
YVO  7 
YV07 
YV07 


YV01  10 
YVO 1  20 

YVO 1  30 

YVO 1  40 

YVO 1  50 

0YV01  60 

YVO 1  70 

YVO 1  80 

YVO 1  90 

YV33 
YV33 
YV3  3 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 
YV33 

YV33 

YV33 

YV33 

YV33 

YV33 

YV18D  10 
YV18D  20 
YV18D  30 
YV18D  40 
YV18D  50 
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941  IF  (CONK  -  0.999999)  942*942.943 

942  CALL  VELLIP  (CONK, EE. EK) 

981  OIF  =  EK  -  EE 

RETURN 

943  IF  (TERM)  944,945,944 

945  DIF  =  1.0 
GO  TO  946 

944  TERM2  =  TERM**2 

ALAM  =  LOGF ( ABSF ( 4 .0*SQRTF ( TERM 2  +  1.01/TERM)) 

OIF  =  ALAM  -  1.0  -  ALAM*TERM2/(4.0*< TERM2  +  1.0)) 

946  EK  =  OIF  +  1.0 
EE  =  -1.0 
RETURN 

END 

FOR 


YV18D  60 
YV18D  70 
YV18D  80 
YV18D  90 
YV18D100 
YV18D110 
YV18D120 
YV18D130 
YV18D140 
YV18D150 
YV18D160 
YV18D170 
YV18D180 
YV18D190 
YV18F  10 


SUBROUTINE  CAMBER  ( DELTA »DN ,H ,N ,NJN ,NZ , P I , TH I CK ,WAKE ♦ WD I FFR , PRAD 
1 ? Jl, THETAB, THETA,  ATTACK, CPRIME,S,Z ) 

DIMENSION  THICKI 1 0 1 > »WD»  FFR (37 ) . PRAD 1 37 )» THETAB ( lO l ) , THETA (37) ♦ 
lCPRIME(lOl ) »CINTl ( 1C1 ) ,CINT2( lOl ) »CINT3( id ) »  VDIFR(lOl)  , 

2S ( 21 ) ,Z ( 21 ) 

N 1  =  N— 1 
RADIAL  =  0.0 
DO  50  1=2. N1 

CALL  INTRL1  ( C INT 1 ( I ) , DELTA, H,  MJN »NZ ,P I , WAKE 

l.WDI FFR, THETA »N»THETAB»VDIFR, I ) 

CALL  MUSK  ( DELTA » ThETa »WD I FFR *NZ ,N ,H»NJN , C I NT3 ( I ) ,TH£TaB* 

1  I , WAKE, PI ,VDIFR) 

CALL  INTRL2  ( N 1 , I , THETAB ,H» 2 ,S t BELT A, Th I CK »N JN » N , C I NT2 ( I ) »P I ) 

GO  TO  (48,50) « J1 

48  CALL  DISCOT  ( THfTAB ( I )♦ THETAB « I ) .THETA , PRAD » PRAD* -120 »NZ , 0 ,RADI AL 

50  CPR I  ME  (  7  -1 )  =  (CINTKI)  +  CINT2  (  I  )  +  C I NT3  (  I  )  -  RADIAL/WAKE  )* 
1SINF(THETA8( I ) )/2.0 

DO  49  1=2. N1 

49  THETAB(I-l)  =  THETAB(I) 

N2  =  N-2 

CALL  DISCOT  (  0 .0  ,C  .0 ,  THETAB  »CPR  I  ME ,  CPRIME  ,-020  »N2  i>0  ?  F I RST  ) 

DO  51  1=1, N2 
NMINUS  =  N-I 
NLESC  =  Nl-I 

THETAB (NMINUS)  =  THETAB ( NLESS ) 

51  CPR I  ME ( NMINUS )  =  CPRIME ( NLESS ) 

THETAB( 1 )  =  O.Q 

CPR I  ME ( 1 )  =  FIRST 

CALL  DI  SCOT  ( P I ,P  I  .THETAB, CPRIME .CPRIME ,-020 »N1 ,0 »CPR IME ( N ) ) 

TANA  =  CPRIME(l)  +  4.0*CPRIME (Nl )  +  CPRIME(N) 

DO  52  I  =2»NJN»2 

52  TANA  =  TANA  +  4.0*CPR IME ( I )  +  2 .0*CPR IME ( I +1 ) 

TANA  =  DELTA*TANA/3.0 

ATTACK  =  57.295780*ATANF(TANA) 

DO  53  I  =  1 » N 

53  CPRIME ( I )  =-CPRIME(I)  +  SINF(THETAB( I ) )*TANA/2»0 
RETURN 

END 

FOR 

SUBROUTINE  INTRL2(N1,I ,THETAB,H,Z,S*DELTA,THICK,NJN,N,SUMI2,PI ) 
DIMENSION  THETAB ( 10 i 1 «F2 ( 101 i .THICK! 101 i ,5THETAi 101 ) ,5(2 1) ,Zi 2ii 
IF{ 1-2 ) 290,290,291 
290  DO  292  J=2»N1 

ZTHETA=.5»(  l.C  t-COSFl  THETAB ( J)  )  ) 


YV320010 

YV320020 

YV32OO4O 

YV32OO5O 

YV320080 

YV320090 

YV320100 

YV320110 

YV320120 

YV320140 

YV320170 
)  YV32O18O 
YV320190 
YV320200 
YV320210 
YV320220 
YV320230 
YV320240 
YV320250 
YV320260 
YV320270 
YV320280 
YV320290 
YV320300 
YV320310 
YV320320 
YV320330 
YV320340 
YV320350 

YV320370 

YV320380 

YV320390 

YV320400 

YV320410 

YV33 


YV33 

YV33 

YV33 
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292  CALL  DISCO T(ZT HE TA.ZTHETA»Z»S»S*-130»21 *0 » STHE TA  ( J ) )  YV33 

291  DO  200  J=2  *N1 

I P ( I  — J  >210*211 *210 

211  F2(J)-0.0  YV33 

GO  TO  200  YY33 

210  PART=C05F (THE TAB ( i ) )-COSF(THETAB( J) )  YV33 

PARTH=.5*H»PART  YV33 

C0NK2=SQRTF( 1 • 0 / < PARTH**2+1 .0  ) )  YV33 

CALL  CWORKI C0NK2»PARTH*EK»EE tDIF )  YV33 

F2( J)=STHETA( J)*CCNK2**3«PART*DIF*SINF(THETAB( J) )  YV33 

200  CONTINUE  YV33 

SUMI2=4.0*F2(N-1)  YV33 

DO  201  J=2  *NJN  *2  YV33 

201  SUMI2=SUMI2+4.0*F2(J)+2.0*F2f J+l ) 

SUM I 2=«25*H**2*DELTA*SUMl2/3»0  YV33 

DO  202  J=  1  *N 

IF(J“1)250*250*251  YV33 

250  F2(J)  =  -STHETAI I ) *S I NF ( THETAB  ( I ) )/( COSF ( THETAB (  I  ) )  -  COSF ( THETABt 

1 J )  )  )  YV33 

GO  TO  202  YV33 

251  IF(N-J>250, 250*252  YV33 

252  I F C I -J 1253*254*253  YV33 

254  F2 ( J )  =  THICK!  I )  +  STHETAU ) «COSF ( THETAB (I ) ) /SI NF { THETABt I)) 

GO  TO  202  YV33 

253  PART=C0SF(THETAB(  m-COSF(THETAB(J) )  YV33 

PARTH=«5*H#PART 

C0NK2=SORTF(l«0/(PARTH**2+1.0) >  YV33 

CALL  EELLIP(C0NK2«EE)  YV33 

F2< J)=(STHETA( J)*S INF (THETABt J) ) *CONK2**3*EE 
1  -  STHETAt I )*SINF(THETAB( 1 ) ) ) /PART 

202  CONTINUE  YV33 

ST0RE=F2(1)+4.0*F2(N-1)+F2(NJ  YV33 

DO  203  J=2  #NJN  1 2  YV33 

203  STORE=STORE+4.0*F2(U)+2.0*F2U+1 )  YV33 

STORE=  DELTA*STORE/3.0  YV33 

SUMI 2=-H* ( SUM  1 2+STORE  J /P I  YV33 

RETURN 

END  YV33 

FOR  YV30 

SUBROUTINE  INTRL1  ( CINTl *DELTA*H »  NJN.NZ.PI *WAKE»WDIFFRYV30 

1>THETA,  N » THETAB*  VDIFR.I)  YV30 

DIMENSION  WDI FFR ( 37 )» THETA ( 37 )» THETAB (101) »F1 ( 101 ) »VDIFR( 101 )  YV30 

IF  (1-2)  202*202*203  YV30 

202  DO  212  J=1»N  YV30 

212  CALL  DISCOT  ( THET«B ( J ) »THET A8 ( J ) *THETA*WDI FFR.WdI FFR  *-020  »NZ  »0  YV3O 

1 »VDI FR ( J ) )  YV30 

203  DO  200  J=1 »N  YV30 

IF  (I-J)  232,231*232 
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231  F1(J)  *  0.0 
GO  TO  200 

232  PART  *  COSF(THETAB< I ) )  -  COSF { THETABt J) ) 
PARTK  =  • 5*PART*H 

CONK1  =  SQRTFt 1.0/(PARTH**2  +  1.0)) 

CALL  CWORK  (CONK1 .PARTH* EK »EE *DI F ) 

F 1 ( J )  =  VDIFR ( J)*PART*CONKi*DIF 
200  CONTINUE 

CINTl  =  Fl(l)  +  4.0*F1(N-1)  +  F1(N) 

DO  250  I -2 *NJN *2 

250  CINTl  =  CINTl  +  4*0*F1(I)  +  2.0*FH1+1) 
CINTl  =-H**2*DELTA»ClNTl/( 12.0*PI*WAKE) 


YV30  120 
YV30  130 


YV30  150 
YV30  160 
YV30  170 
YV30  180 
YV30  230 
YV30  240 
YV3O  250 
YV30  260 
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RETURN 

END 

FOR 

SUBROUTINE  INTGRL  (ANGLE  .CONSK » COSl , EE » FB , F IRST » 
1NNN,  THETA1 ♦  THETAS ,JJJ ,M » H ) 

DIMENS (ON  CONSk { 1 0 1 ) »THETA5 ( lO l ) *EE( 1 01 ) * F 3 ( lOl ) 
LOOP  FOR  CALCULATION  OF  INTEGRAND  VALUES  FOR  INTEGRAT 

RESPECT  TO  THETA  5 
DO  1-Oi  J= J  JJ »NN » JJ J 
IF  ( J-JJJ)  800  *800,810 
800  THETAS(J)  =  FIRST 
THETA5  VARIES  FROM  0.0  TO  PI 
GO  TO  811 

810  THETA5 ( J)  =  THETAS(J-l)  +  ANGLE 

811  DEN  =  COSF ( THFTA5 ( J ) )  -  COSl 

IF  (ABSF(DEN)  -  0.000001)  850*850*856 
850  FB ( J )  =  0.0 
GO  TO  101 

856  CONSK(J)  =  SQRTF(l.0/(0.25*(H*DEN)**2  +  1.0)) 
CALCULATION  OF  ELLIPTIC  CONSTANT  K. 

CALL  EELLIP  ( CONSk ( J )* EE ( J ) ) 

FB ( J  )  =  (1.0  -  1  .C /CONSK ( J )  *  EE(J) )  * ( S I NF ( THETA5 
PARTIAL  INTEGRAND  VALUE  FOR  INTEGRATION  WITH  RESPECT 
101  CONTINUE 
RETURN 
END 
FOR 

SUBROUTINE  MUSK  (DELTA*  THETA *WD1 FFR »NZ *N *H 

1  *  I .WAKE, PI ,VDIFR) 

DIMENSION  THETA (37) ,WD I FFR ( 37 J  » F3 (10 1 ) * VD I FR ( 10 1 ) 
DO  10  J=1*N 
IF  (I-J)  Il9l2.ll 
12  ABOVE  =  THETAB(I)  +  .01 

CALL  DISCOT  ( ABOVE , ABOVE ♦ THETA ,WD I FFR ,WD I FFR ,-020 
BELOW  =  THETAB ( I )  -  .01 

CALL  DISCOT  ( BELOW .BELOW , THETA ,WD I FFR ,WD I FFR ,-020 
F  3 ( J )  *  (TAB  -  FBEL)/( .02*S INF ( THETAB ( I ) > > 

GO  TO  10 

11  DENOM  =  COSF ( THETAB ( I))-  COSF (1 HET AB ( J ) ) 

CCONK  =  SQRTF( 1 .0/ ( 0. 25* ( H*DENOM ) **2  +  1.0)) 

CALL  EELLIP  (CCONK, EE) 

F3 ( J )  =  (VDIFR(J)*CCONK*EE  -  VD I FR ( 1) ) /DENOM 
10  CONTINUE 

CINT3  =  F3(l)  +  4.0*F3(N-1)  +  F3(N) 

DO  40  J=  2  ,N JN  *  2 

40  CINT3  =  CINT3  +  4.0*F3(J)  +  2.0*F3(J+1) 

CINT3  =  DELTA*CINT3/(6.0*PI*WAKE) 

RETURN 

END 
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( J) ) ) **2/DEN 
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,NZ ,0 , FAB ) 


,NZ  ,0 ,FBEL ) 
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